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Summary 

A new quasi-one-dimensional Godunov code for mod- 
elling two-stage light gas guns is described. The code is 
third-order accurate in space and second-order accurate in 
time. A very accurate Riemann solver is used. Friction 
and heat transfer to the tube wall for gases and dense 
media are modelled and a simple nonequilibrium turbu- 
lence model is used for gas flows. The code also models 
gunpowder burn in the first-stage breech. Realistic equa- 
tions of state (EOS) arc used for all media. The code was 
validated against exact solutions of Riemann’s shock-tube 
problem, impact of dense media slabs at velocities up to 
20 km/sec, flow through a supersonic convergent- 
divergent nozzle and burning of gunpowder in a closed 
bomb. Excellent validation results were obtained. The 
code was then used to predict the performance of two 
light gas guns (1.5 in. and 0.28 in.) in service at the Ames 
Research Center. The code predictions were compared 
with measured pressure histories in the powder chamber 
and pump lube and with measured piston and projectile 
velocities. Very good agreement between computational 
fluid dynamics (CTO) predictions and measurements was 
obtained. Actual powder-burn rates in the gun were found 
to be considerably higher (60-90 percent) than predicted 
by the manufacturer and the behavior of the piston upon 
yielding appears to differ greatly from that suggested by 
low-strain rate tests. 

1 Introduction 

In 1993, after being shut down for a number of years, the 
radiation ballistic range facility at the NASA Ames 
Research Center was reactivated. Tests were performed in 
1 993—94 using Ames' 0.28 in. and 1.5 in. caliber two- 
stage light gas guns. These tests studied, primarily, dam- 
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age caused by simulated space debris impacts on space 
station wall segments and space shuttle tiles. Ames has 
long served as an agency source of expertise in light gas 
gun launcher research, development and testing (ref. I). 
For the space debris impact tests, an increase in the gun 
muzzle velocity was desired to allow a more complete 
simulation of the possible velocity range of space debris 
impacts. Such a muzzle velocity increase must be accom- 
plished with great care to avoid excessive gun erosion and 
the overstressing of the gun or the launch package. It 
would also be desirable to reduce the maximum gun and 
projectile base pressures and gun erosion, while maintain- 
ing muzzle velocity. Experimental gun development is 
very costly in time and money (there can be eight or more 
gun operating parameters which can be varied) and can 
carry considerable risk of major damage to the facility. A 
well-validated, user-friendly computational fluid dynam- 
ics (CFD) gun code capable of guiding the selection of 
gun operating parameters to safely and economically 
achieve increases in muzzle velocity and/or reductions in 
maximum gun pressures and gun erosion is needed. The 
code described herein has been used to help increase the 
muzzle velocity of the Ames 1 .5 in. gun (for a fairly 
heavy projectile) from 6.6 to 7. 1 km/sec, while, at the 
same time, reducing gun erosion by 50 percent. 

Figure 1 shows a sketch (not to scale) of a representative 
two-stage light gas gun. From left to right, we see the 
powder chamber, the initial position of the plastic piston 
in the pump tube, the very long pump tube, initially Filled 
with hydrogen or helium, the contraction section, the 
diaphragm just behind the projectile, the projectile in its 
initial position and the gun barrel (or launch tube). Upon 
ignition of the powder charge, the piston is accelerated to 
a velocity of the order of 800 m/sec. The piston travels 
down the pump tube, greatly compressing and healing the 
working gas. The working gas can be volumetrically 
compressed by a factor of the order of 1000, reaching 
pressures of roughly 10,000 bar and (theoretical) tempera- 
tures as high as 3000 K. At certain point in the 
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compression process, the diaphragm ruptures and the 
projectile begins to accelerate down the gun barrel. 

The great advantage of the two-stage gun over a single- 
stage powder or gas gun is as follows. For the single-stage 
gun, it is difficult to obtain a sound speed of the driving 
gas which is much greater than 1—1.4 km/sec. This usually 
limits the muzzle velocity of such a gun to a maximum of 
2.5-3 km/sec. With the two-stage gun technique, the 
energy of the powder gas is transferred, with a reasonably 
high efficiency, to a much smaller mass of low molecular 
weight working gas. The sound speed of the working gas, 
at maximum compression, can be as high as 4^4.5 km/sec. 
This permits the two-stage gun to reach muzzle velocities 
of 7-8 km/sec fairly routinely. With very light projectiles, 
and driving the gun hard, velocities as high as 
10-11 km/sec can be obtained. 

All media in the gun (powder/powder gas, piston, working 
gas, and projectile) are subject to very high pressures and 
temperatures and therefore should be modelled with 
equations of state (EOS) which arc accurate under these 
conditions. As an example, for the gases, the ideal gas 
EOS is usually very inaccurate in the setting of two-stage 
guns. The piston and projectile may show high strain rate 
behavior which is very different from that which would he 
expected based on available low-strain rate data. Friction 
and heat transfer of gases and solids to the tube wall arc 
very important, can cause reductions (refs. 2 and 3) in 
muzzle velocity of 10-20 percent, and should be carefully 
modelled. It would be desirable to assess to what extent 
nonequilibrium turbulence in the gas flows may affect the 
skin friction and heat transfer to the tube walls. 

The two-stage light gas gun was not the only quasi-one- 
dimensional How problem which needed to be modelled 
at Ames at this point in time. There was also interest in 
high-explosive detonations, high explosive driven ram 
accelerators and hypcrvelocity impact. These phenomena 


require accurate modelling at pressures from 300 kbar to 
10 Mbar and higher pressures, well beyond the 
representative maximum pressures in light gas guns 
(10-15 kbar). It would be desirable to obtain a code which 
was capable of dealing with these very-high-pressure 
problems as well as for modelling light gas guns. 

A number of CFD codes (refs. 3-9) for one- and two- 
stage guns were examined. The IBHVG2 code (ref. 5) is 
limited to modelling single-stage guns. The code by 
Charters and Sangster (ref. 3) is a two-stage gun code, but 
does not include friction and heat transfer and uses EOS 
which have limited ranges of applicability. It is known 
(refs, 2 and 3) that friction and heat transfer effects can 
cause muzzle velocity reductions of 10-20 percent. The 
ALE code (ref. 9) docs not calculate the entire gun cycle, 
but rather starts with a given piston velocity history and a 
point mass piston. Powder-burn and piston deformation 
are not treated. The gas friction model does not include 
Mach number and wall-temperature effects. For the EOS 
of hydrogen, a 28-coefficient fit to tabulated data is used. 
The MOOREA code (refs. 7 and 8) uses a van der Waals 
volumetric EOS for the pump tube gas and powder gases 
and a segmented thermal EOS for the pump-tube gas. The 
latter EOS is composed of three regions, each with C v 
varying as AjT + Bj, where C v is the specific heat at con- 
stant volume, T is temperature and Aj and Bj arc con- 
stants, different for each region. The friction model for the 
piston and projectile either uses a constant force at the 
base of the piston or projectile or increases the mass of the 
piston or projectile a few percent. The LLG3 model 
(ref. 4) uses a volumetrically and thermally perfect EOS 
for the pump tube gas. The gas friction model does not 
include a wall-temperature effect and wall heating from 
piston and projectile friction do not appear to be included. 
The HVML89 model (ref. 6), developed from the earlier 
model of reference 2, uses the van der Waals EOS for the 


2 




pump-tube and powder gases and an isothermal EOS for 
the piston material. 

None of these codes included nonequilibrium gas turbu- 
lence effects and the EOS limitations would preclude 
them from producing accurate results at megabar pres- 
sures. In addition, all of the code examined failed to 
include one or more additional phenomena which it was 
desired to model. The authors emphasize that there is 
nothing wrong with the codes discussed above when used 
for the purposes for which they were written and they 
have, in fact, yielded a large number of very useful 
results. Rather, they do not easily lend themselves to 
model all the phenomena we wished to study in two-stage 
light gas guns and also to study megabar pressures (at 
least, not without substantial modifications). Making sub- 
stantial changes deep within another worker’s program 
can be an extremely lengthy and, at times, cost-ineffective 
process (although one may learn to run a well- 
documented code written by another worker rather easily). 
For these reasons, taken together, it was decided to 
construct a new code. 

The new code presented here uses the Godunov (ref. 10) 
technique. Briefly, the Godunov technique solves a 
Riemann problem at every cell boundary at every timestep 
in order to calculate the flux between the two cells in 
question at the lime in question. Figure 2 shows an 


x-l wave diagram of such a Riemann problem. The 
abscissa is distance (perpendicular to the cell boundary) 
from one cell to the other and the ordinate is time starting 
from the beginning of the timestep. The initial conditions 
on the two sides of the cell boundary at the beginning of 
the timestep are conditions I and 4. The dynamics of the 
Riemann problem gives rise to a pressure wave travelling 
into the medium of the left cell, an interface between the 
two media and a pressure wave travelling into the medium 
of the right cell. These waves and the interface are shown 
in figure 2. The conditions in the two new zones 2 and 3 
can be calculated to many different degrees of accuracy, 
as required. The flux is then calculated according to the 
zone in which the cell-boundary resides during the 
timestep in question. For the case of figure 2, if the cell 
boundary were motionless, its path in the x-t diagram 
would be vertical and it would reside in zone 2. The 
fluxes would then be evaluated simply and directly from 
the media conditions of zone 2. The technique is exceed- 
ingly flexible. No assumptions about the EOS of the 
media are necessary. It is not necessary to explicitly insert 
any artificial viscosity. Extrapolations (or interpolations) 
to the cell boundary in question can be first, second, third, 
or even higher order accuracy. The Riemann solver can be 
first, second, or higher order accurate and its order of 
accuracy can readily be changed to suit the problem in 
question. 



Figure 2. x-t wave diagram for Riemann problem between two cells . Initial conditions at beginning of timestep are those in 
regions 1 and 4. Waves progress into both these regions and an interface separates the media initially in regions 1 and 4. 
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The features of the new higher-order Godunov code pre- 
sented herein include the following: 

( 1 ) It models the complete two-stage gun cycle, from 
the start of powder burn, to the projectile exiting the 
muzzle. 

(2) It uses the Godunov technique, outlined above, 
which has been shown (ref. 10) to be very robust for 
dense media and extremely high velocities. It is third 
order accurate in space and second-order accurate in time. 
The Riemann solver used is of very high accuracy, being 
exact for shocks and using a very accurate power law 
integration scheme for expansion waves. 

(3) It uses realistic EOS for all media. The EOS for 
solids includes the effects of both density and energy on 
pressure and would remain valid for megabar pressures. 
Simple models for high strain rate effects can be intro- 
duced into the modelling of the piston and the projectile. 

(4) Friction and heat transfer from gases and dense 
media to the tube w r all are included. The gas-phase 
friction and heat-transfer model includes compressibility 
and wall-temperature effects. A simple model for 
nonequilibrium gas turbulence is included. 

(5) The basic algorithms are capable of producing 
accurate results at pressures well into the megabar range. 
Thus, the code could be used to model high explosive 
detonations or hypervelocity impact phenomena. During 
the validation of the present code using exact analytical 
solutions (see sec. 4.2), accurate CFD solutions were 
obtained at pressures up to 8 megabars. A earlier version 
of the code (ref. 10) has produced accurate solutions at 
pressures up to 10,000 megabars. 

(6) The code is well-documented and user-friendly 
versions exist for Cray and HP computers. The documen- 
tation includes descriptions of all the physics modelled, 
the algorithms used and approximations made. A manual 
is available, giving complete instructions for the use of the 
code. 

This technical memorandum will discuss all aspects of the 
higher-order Godunov code. The governing equations, 
source terms and boundary conditions (be) will be dis- 
cussed along with the accompanying numerical methods. 
The code will be validated, in part, by comparison with a 
number of exact analytical solutions. Further validation 
will be shown by comparison with recent, extensive 
experimental data from Ames’ 0.28 in. and 1.5 in. caliber 
light gas guns. The code has been found to be very useful 
in a theoretical and experimental study (ref. 1 1) of the 
technique of inserting an extra diaphragm into the pump 
tube of a two-slagc gun in order to increase gun perfor- 
mance. The code has been used to guide modifications of 


the operating conditions of the Ames 1.5 in. gun to yield 
muzzle velocity increases (fora moderately heavy projec- 
tile) from 6.6 to 7.1 km/sec and, at the same time, a 
50 percent decrease in gun erosion (ref. 12). Further, it 
has yielded some very interesting insights into the burning 
of gunpowder in the first-stage breech and the apparently 
anomalous yielding behavior of the pump tube piston. 
These last two results will be presented herein. 

Support for DWB by NASA (Contract NAS-2- 1 403 1 ) to 
Eloret is gratefully acknowledged. 


2 Formulation 

In this section, wc present the formulation of the problem. 
The governing gasdynamic equations are presented in 
section 2. 1 . The various options available for the EOS of 
the media are presented in section 2.2. Sections 2. 3-2. 6 
present the formulation of the various source terms which 
appear in the equations. These include the wall-pressure 
term, friction and heat-transfer terms for gases and dense 
media and the gunpowder combustion terms. 


2.1. Governing Equations 

The quasi-one-dimensiona! gasdynamic equations used, 
written in conservation form, arc 


A 


3U 

3t 


3(FA) 3A 

~^r +H ih 


+ AW + tcDV 


(!) 


where the state vector is given by 


U = (p,pu,e t ,pm|,pm 2 ) (2) 


the flux vector is given by 


F = (pu,pu 2 


+ p,u(e t + p),pum] ,pum2,...) 


(3) 


the wall-pressure component of the source term is given 
by 


H = (0,p, 0,0,0,...) 


(4) 


the chemical-reaction component of the source term is 
given by 
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W = (0,0,0, W|,W2,...) (5) 

and the wall-friction and heat-transfer component of the 
source term is given by 

V = (0,T,q,0,0,...) (6) 

where A denotes cross-sectional area, D, tube diameter, p, 
density, u, velocity, e, the internal energy per unit mass, p, 
pressure, mj, mass fraction of the ith component, 
q = p(e + u 2 /2), wj, the species mass generation rate per 
unit volume for the ith component, -x, the wall-shear 
stress and -q, the wall-heat flux. The speed of sound, c, 
temperature and pressure are obtained from the EOS 
described below. The code uses a finite volume formula- 
tion with state variables calculated at the center of each 
computational cell. 


experimental data of reference 10 up to megabar pres- 
sures, and 

(3) its inclusion significantly complicates the solu- 
tion of the EOS in the code and slows down the code 
appreciably. 

From the experimental shock Hugoniot data of ref- 
erence 16 for a given medium, a good determination of 
the constants in this EOS can be made. Option #4 is a 
two-dimensional tabulated EOS of the form T(p,e), p(p,e). 
From the tabulated grid-point values, bi-cubic interpola- 
tions are made in logarithmic space to obtain values of T, 
p, Pp,e> ar| d p c ,p at tbe desired point. The sound speed is 
then calculated from p, p, pp <c , and p e p. For the hydrogen 
gas used in the current CFD code, the EOS was con- 
structed as follows. First, equilibrium calculations were 
made for point hydrogen molecules (and atoms) over the 
full required pressure and temperature ranges. Then, a 
molecular volume term was added which follows the 
Zefdovich and Raizer (ref. 15) cold pressure-volume 
relation for dense media. 


2.2 EOS 

The code offers a number of EOS options. We will 
describe them below. (The label numbers given refer to 
those used in our CFD code.) Option #1 uses the Abel 
(ref. 13) volumetric EOS, 


p(v-b) = ^ 
m 


(7) 


where p is pressure, v, specific volume, b, molecular vol- 
ume, R u , the universal gas constant, T, temperature and 
m, molecular weight. The enthalpy is determined assum- 
ing constant specific heat. This option can be used with 
single or multiple species. Option number #2 also uses 
equation (7), but determines the enthalpy with variable 
specific heats using data from the Joint Army-Navy- 
NASA-Air Force (JANNAF) tables (ref. 14). Again, this 
option can be used with single or multiple species. Option 
#3 is Zefdovich and Raizer’s three term dense media 
EOS (ref. 15) with the third term neglected. This term 
was dropped because: 

(1) it is an electron temperature term which is very 
small for the dense media conditions of interest here, 

(2) calculation of shock Hugoniots for the piston and 
projectile material using two terms only of the Zefdovich 
and Raizer EOS shows excellent agreement with the 


P ~ A 


Po y 


-1 


( 8 ) 


where p 0 is a reference density and A and n arc empirical 
constants for the medium in question. The constants in 
equation (8) for hydrogen were determined by fitting the 
experimental shock Hugoniot data (ref. 16) for liquid 
hydrogen to the Zefdovich and Raizer’s three term dense 
media EOS (ref. 15) with the third term neglected. It can 
be shown that this yields an EOS for hydrogen which 
compares rather well with the tabulated SESAME 
(ref. 17) EOS. There arc about twice as many grid points 
per decade in our gridding than for the SESAME EOS, 
which was designed for an extremely wide range ol densi- 
ties and temperatures, i.e., densities from 5 x 10~^ to 
1750 gm/cnv 1 (for hydrogen) and temperatures from 0.04 
to 30,000 electron volts. The higher density of grid points 
considerably improves the accuracy and smoothness of 
interpolated results. 

Two options are available for two-phase EOSs for 
the gunpowder/powder gas regions. Option #7 uses 
the tabulated EOS for the gas phase and a simple model 
with constant density and internal energy for the solid 
phase. Option #8 uses the option #1 (with one species 
only) for the gas phase and the constant density, constant 
internal energy model for the solid phase. In the gun 
solutions presented herein, option #8 was used lor the 
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gunpowder/powder gas regions, option #4 for the piston 
and projectile and option #6 for the hydrogen gas. 

2.3 Source Terms-Wall-Pressure Term 

When there is an area change over the length of a compu- 
tational cell, a momentum source term is required in the 
equations. A simple way to evaluate this term would be as 
follows. 


F w = £i y El -AA (9) 


where F w is the wall -pressure force (source) term, AA is 
the area change and pj and p r are the pressures at the left 
and right ends of the cell. The average pressure acting 
over the wall is assumed to be (p| + p r )/2. A more sophis- 
ticated technique can be used to evaluate F w , as follows. 
Parabolic fits arc made to the channel area profile and the 
pressure variation along the channel. These can then be 
integrated over the cell length to provide a more accurate 
estimate of F w . In code proof tests in a supersonic con- 
vergent-divergent nozzle, the more sophisticated evalua- 
tion ofF w yielded significantly lower errors and, hence, 
has become a permanent part of the code. 

2.4 Source Terms-Gas-Phase Friction and Heat- 
Transfer Terms 

The calculation of these terms is presented in some detail 
in Appendix A and, hence, will only briefly be outlined 
here. We start with the well known (refs. 18 and 19) skin 
friction correlations for pipe flow. We base a correction 
for the effects of Mach number and the difference 
between the wall temperature and the average stream 
temperature on the correlations developed by van Driest 
(ref. 20). We have developed a "reference temperature” 
technique which fairly accurately reproduces the varia- 
tions of skin friction coefficient, f, with Mach number and 
the (wall tempcrature)/(strcam temperature) ratio given by 
van Driest’s correlations. The turbulent component of f is 
then corrected for nonequilibrium turbulence effects as 
described in Appendix B and briefly, in the next para- 
graph. With f determined, the wall friction momentum 
term follows directly. We then use a form of Reynolds' 
analogy (ref. 21 ) to estimate the wall-heat transfer energy 
term. For regions with gas and solids (gunpowder/powder 
gas), the wail friction and heat-transfer terms arc based on 
the gas properties only. This two-phase flow condition 
applies only in the upstream part of the pump lube for the 
first part of the solution before all the powder is burned. 
Currently, T w is prescribed and held fixed at 300 K. At a 


later time, the code will be modified to allow for wall- 
heating effects. 

A simple model (App. B) was developed which assumes 
that the nonequilibrium turbulence kinetic energy (TKE) 
relaxes towards the equilibrium value for the flow in ques- 
tion with an e-folding length which is a certain number of 
tube diameters. (The c-folding length is the distance along 
the tube over which the difference between the nonequi- 
librium and equilibrium TKE would drop to 1/cth of its 
original value in a steady flow with constant cross- 
section.) The range of Reynolds numbers (Re) for hydro- 
gen flow in the pump lube is typically 3 x 10^ to 3 x 10 7 . 
Detailed turbulence measurements (refs. 22 and 23) for 
fully developed pipe flow were found at a maximum Re 
of 5 x 105 This Re is within our range, but towards the 
low end of it. However, turbulent pipe flow does not 
appear to change very rapidly with Re over the Re range 
of interest (at least over the range 3 x 10^ to 3 x 10^ 
reported in ref. 24). Hence, a rough estimate of the 
number of pipe diameters required for relaxation of the 
turbulent kinetic energy towards the equilibrium value 
was obtained from the data of references 22 and 23. This 
number (3.27 diameters) was used in our model. By writ- 
ing a simple TKE equation with this relaxation term and 
the appropriate convection terms across cell boundaries, 
the TKE in any cell at any time can be calculated. With 
values for the nonequilibrium and equilibrium TKEs in 
hand, the turbulent portion of the skin friction coefficient 
is multiplied by a term (TKE none q U j|/TKE C q U ji)° ^ to 
reflect the increased transport of momentum and energy 
due to nonequilibrium turbulence. Further details arc 
given in Appendix B. 

2.5 Source Terms-Solid-Phase Friction and Heat- 
Transfer Terms 

The calculation of these terms is also presented in some 
detail in Appendix C and will only briefly be outlined 
here. Our technique is a variation of that presented in ref- 
erence 4. First, the normal stress on the tube walls is cal- 
culated from the axial pressure which appears in the CFD 
code using a simple elastic-plastic model. This stress is 
first calculated, assuming clastic behavior, taking into 
account (1) the initial jamming of the solid into the tube 
and (2) the axial (CFD) pressure term. For the piston, the 
initial jamming is produced by cooling the piston in a 
freezer (which causes it to contract), inserting it into the 
pump tube and allowing it to return to room temperature. 
As the piston temperature rises, it expands and generates 
jamming stresses. If the calculated difference between the 
axial and normal stresses is found to be greater than the 
yield stress, this difference is simply set to the yield stress. 
This latter procedure is the plastic part of the modelling. 
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The friction coefficient between the solid and the walls, as 
a function of velocity, is modelled from the experimental 
data of reference 25, following the techniques of refer- 
ences 26-28. The friction force on the wall is then calcu- 
lated and limited by the estimated shear yield stress of the 
solid, following reference 29. The momentum source term 
follows directly. The energy source term is simply taken 
to be 1/2 x (friction source term) x (velocity). This term is 
a loss of energy from the solid to the wall due to frictional 
heating. The other half of the heat generated by the fric- 
tional work is assumed to flow to the solid and therefore 
docs not represent an energy loss from the solid. 

2.6 Source Terms-G unpowder Burn 

The powder grains are of the standard form of circular 
cylinders with a number (usually 0, 1 , or 7) of circular 
perforations parallel to the outer cylinder axis. The linear 
surface regression rate, r, of the powder is taken to follow 
the usual ballistics expression (sec, for example, ref. 30) 

r = ap n (10) 


the specific energy release of the powder upon combus- 
tion known, the necessary powder burn source terms for 
energy deposition and change of powder into powder gas 
can readily be evaluated. The amount of powder burned 
per cell per timestep is calculated using a predictor- 
corrector method, so it is second-order accurate in time. 

The unburned powder and the powder gas are assumed to 
move together — there is no relative motion between the 
two phases and, hence, no erosive burning effects due to 
relative gas-powder velocity differences. The powder is 
also assumed to be distributed at a uniform density 
throughout the first-stage breech. This avoids the produc- 
tion of spurious slosh wave effects. These were deliber- 
ately accepted limitations made to allow the code effort to 
concentrate on two-stage effects, that is, piston, hydrogen, 
and projectile dynamics. The code thus is not suitable for 
the detailed study of single-stage military-type weapons 
where relative motion between the gas and the powder 
and erosive burning effects can be important. From the 
comparison between experimental and theoretical 
powder-chamber and pump-tube pressure histories and 
piston velocities (see secs. 5.3 and 5.4), the above 
assumption will prove to be very good for the present 
purposes. 


where p is the pressure and a and n arc constants given by 
the maker of the powder or by military testing laborato- 
ries. The specific energy released upon burning of the 
gunpowder, AE, is calculated from 


AE = E(powder gas) - E(solid powder) 

(id 

E(powder gas) = 

(12) 

Y-l 


I = Impetus = RuTa 
h M 

(13) 


where E(solid powder) is the internal energy of the solid 
powder, y is the specific-heat ratio of the powder gas, R u 
is the universal gas constant, is the adiabatic flame 
temperature of the powder and M is the mean molecular 
weight of the powder gas. I (the ‘impetus”), y, and M are 
given by the powder manufacturer. E(solid powder) is 
estimated at room temperature from the specific heat of 
the powder. The density and energy of the unburned pow- 
der are assumed to remain constant at their room tempera- 
ture values. Thus, no heat transfer to the powder is 
considered. 

With the initial shape of the powder grains known, the 
linear regression rate calculated using equation (10) and 


3 Numerical Method 

The numerical method will be described in this section. 
The method of updating of the cell-center state variables 
in time is presented in section 3. 1 . The extrapolations, 
interpolations and limiters used to calculate the primitive 
variables at the cell boundaries are described in sect- 
ion 3.2. The Riemann solvers used to calculate the fluxes 
at the cell boundaries are presented in section 3.3. 
Calculations of the conditions at the zone boundaries arc 
described in section 3.4. The multiple zoning techniques 
used, including the selection of the distributions of cell 
sizes in the various zones, are described in section 3.5. 
Finally, the techniques used to model jamming or bounc- 
ing of the pump-tube piston are presented in section 3.6. 

3.1. Code Advancement in Time 

To calculate the fluid dynamic part of the time step, the 
code employs an explicit MacCormack predictor-corrector 
differencing scheme (ref. 31) which is second-order 
accurate in time. For this purpose, the solution to equa- 
tion (I) is advanced in time for the predictor and corrector 
as follows. 


\j n+1 y n + l 


U"V n 



+ F S 


At 


(14) 
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Fl At-FgAg +F"' 


At (16) 


Fg’ = (0, pAA + xA w ,qA w ,0,0 )"' 


(17) 


pn , F n + F n+I 


(18) 


The best performance was found when fer was in the 
range 03-0.5. The results presented herein were calcu- 
lated with f cr equal to 0.5. This safety device is invoked 
only very rarely, when very strong shock waves are 
involved. 

After the above second-order fluid dynamic advancement 
in time is performed, a second-order time accurate 
“chemistry in a closed box” powder-burn calculation is 
made, as described in section 2.6. For the powder burn, 
the solution to equation (1 ) is advanced in time for the 
predictor and corrector as follows. 


U n +l v n =u n V n +F s n b At 


(19) 


where 

U = state vector 

V - cell volume 

Al.Ar = left- and right-cell boundary areas 

AA = Ar - Al 

A w = projection of cell wall area normal to axis 

F l ,Fr = left- and right-cell boundary fluxes 

F$ = fluid dynamic source term 

At = time step 

q = wall heat flux 

T = wall shear stress 

p = “average” cell pressure 

n superscript denotes conditions at start of step 

n + 1 superscript denotes conditions at end of predic- 
tor step 

n+l superscript denotes conditions at end of corrector 

step 

n 1 superscript denotes average of the fluxes evaluated 
at the n and n + l time levels 

For details of the calculations of the pAA, tA w . and 
qA w terms, sec sections 23-2.5 and Appen dice s A-C. 
When the computational grid slides, V n * V n+ * * V n+ * 
and the areas must he carefully averaged values. Other- 
wise, spurious source terms can be generated. If the cell 
center e or p values calculated from equations (14) and 
(16) arc less than a given fraction (f cr ) of the value at the 
beginning of the timestep, they are reset to that fraction. 


F" b = vn (°,°, wA E b ,-w,w) n (20) 

U n+, V n =U n V n +F s n b At (21) 

F s n b = V " (0, 0, wAE b , - w, w) n (22) 

where 

Fsb = P owc ^ cr burn source term 

w = mass of powder consumed per unit volume, 
per time 

AEf> = energy released upon powder burn, per unit 
mass of powder 

In the formulation of equations ( 1 9)— (22), only two 
species are considered, unburned gunpowder and powder 
gas. Also, since the chemistry takes placed in a fixed, 
closed box, the volume, V n , is the same at all time levels; 
it is included in equations ( 1 9)— (22) simply to place them 
in the same form as equations ( 1 4)— ( 1 8). Since these two 
second-order procedures are completely separated in the 
time step, the interaction between them is only first-order 
accurate in time. For the current purposes of gunpowder 
burn only, this does not cause any difficulties. Since the 
method is explicit, the von Neumann stability criterion 
(ref. 32) that the CFL number be less than one must be 
applied. Most of our results were obtained at a CFL num- 
ber of 0.7. 

3.2 Determination of Cell Boundary Values 

The fluxes at the cell boundaries are calculated from val- 
ues of the primitive variables (p, u, c, and mj) on the two 
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sides of the cell boundary. The cell-boundary values are 
obtained from extrapolation or interpolations of the cell- 
center values. All extrapolations and interpolations arc 
made by cell numbers, taking no account of the physical 
sizes of the cells. Most of the calculations of internal cell- 
boundary values are with third order extrapolations and 
interpolations. Six sets of variables with which to extrapo- 
late and interpolate were tested before deciding on (p, u, 
e, and mj). In these six sets of variables, the second vari- 
able was cither u or pu and the third variable was either e, 
pe, or C(. In extensive testing with the Reimann shock 
tube problem at a pressure ratio of 1000, the primitive 
variables gave the best combination of minimum oscilla- 
tion size (about 1-2 percent) and minimum spurious 
entropy production and heating. (Small, rapidly damped 
oscillations of amplitudes of 1-2 percent sometimes 
appear in the solutions in the neighborhood of strong 
shock waves.) 

Figure 3 illustrates the extrapolation and interpolation 
procedures. We consider the right-side boundary value of 
a primitive variable at the boundary 1-2 between cells I 


and 2. The cell-center values arc indicated by the dots for 
cells 1 to 4. At the boundary 1-2, third-order extrapolation 
gives the value 3e, third-order interpolation gives the 
value 3i and second-order extrapolation gives the value 2, 
all indicated by open circles. The choice of boundary val- 
ues used in the present code is as follows. If the 3i value is 
between the 3e and 2 values, it is used as the final 
cxtrapolation/intcrpolation value. If the interpolated value 
is outside the range of the 3e and 2 values, whichever of 
the 3e and 2 values is closest to the 3i value is taken as the 
final extrapolation/interpolation value. So, for the case 
shown in figure 3, the 3c value would be selected. 

For the cell boundary one cell removed from a zone 
boundary, extrapolation/interpolation from the zone 
boundary towards the cell boundary is generally of second 
order only. This is because only one, rather than two, 
ghost cell is used beyond each zone boundary. The reason 
for dropping to second order at this location is that when 
constructing ghost cells across a boundary between two 
very different media, there is a considerable degree of 
freedom even in constructing a single ghost cell. We were 



Cell number 

Figure 3. Wustration of extrapolation and interpolation from cell-center to cell-boundary values. 
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able lo find reasonable, robust solutions for this problem 
with single ghost cells. However, for two ghost cells, the 
degree of freedom becomes much more severe, and it was 
judged better not to broach the problem, but rather to 
remain with a single ghost cell. 

In addition, for some simple inlet and outlet conditions, 
such as constant condition inlet or constant pressure out- 
let, simple first-order be were used at the inlet or outlet 
zone boundaries. These first-order conditions were used 
only for proofing of the code and were not used for any 
gun solutions presented herein. All of the gun solutions 
presented herein used third-order extrapolations/ 
interpolations, except as described in the preceding para- 
graph. The code has, however, options to be run with 
second- or first-order extrapolations and has occasionally 
been so run, when problems occurred during the initial 
code development and debugging. 

After the cxtrapolations/interpolations of the primitive 
variables have been done to the cell boundaries, they are 


limited to prevent them from being outside the range of 
the cell-center values straddling the cell boundary in ques- 
tion. This procedure is illustrated in figure 4. The abscissa 
is the original extrapolated/interpolated value and the 
ordinate is the final limited value. The diagonal line 
through the origin is at a 1:1 slope. The solid line repre- 
sents the limited primitive variable. Basically, the limited 
value is taken to be the extrapolated/interpolated value if 
it is between the cell-center values, but if it is outside the 
range of the cell-center values, the nearest cell-center 
value is substituted instead. A refinement has been added, 
however, as shown in figure 4. Instead of using abrupt 
breaks at the change points between the use of full 
extrapolated/interpolated values and the use of cell-center 
values, we have found superior code results to be obtained 
if one uses parabolic blending in these regions, as shown 
in figure 4. The size of the parabolic blending regions is 
much exaggerated in figure 4 (by about a factor of 6) for 
clarity. A value of 5/L of 0.025 was found to work well. 



Figure 4. Illustration of parabolic blending of the limiting of extrapolated/interpolated cell-boundary values. 
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The above limiting processes were found to be insuf- 
ficient for calculations of very strong waves. A "strong 
wave limiting” (SWL) process, described below, was 
found to be necessary. Consider the extrapolation/ 
interpolation procedures from the right hand direction to 
obtain the cell-boundary conditions between cells 1 and 2 
in figure 1. Across each of the cell boundaries involved in 
the extrapolation/interpolation procedures, i.e., bound- 
aries 1-2, 2-3, and 3-4, we assess the strength of the 
waves present by calculating the following three 
parameters, pi, 02, and P3. We then find the largest 


Pi=- 


|U| -u r | 


min 


(C|-C r ) 


(23) 


P2 



|Pl - Pr 

min| 

2 2\ 
[P|C| ’Pr c r j 


(24) 


P3 = 


|p| c l ~ Pr c rl 
min(p|Ci,p r c r ) 


(25) 


of the nine P values (three types of P values, each at three 
different cell-boundary locations) and denote it by Pmax- 
(In eqs. (23)-(25), c denotes the sound speed and the sub- 
scripts 1 and r denote conditions on the left- and right- 
hand sides of the boundary.) We use this P max value to 
decide if strong waves are present in the extrapolation/ 
interpolation zone and hence, whether further limiting 
towards first-order values is required. If Pmax * s greater 
than a critical value, otj, the extrapolated value is replaced 
the corresponding first-order cell-center value. If Pmax 
less than a second critical value, 0 C 2 , the higher-order 
cxlrapolatcd/interpolated value is used without further 
modification. If p max is between ot| and ct 2 , a linear blend 
of first- and higher-order values is used. This linear blend 
assures that there is always a smooth, progressive change 
from higher order to first oeder as the wave becomes 
stronger. The values for a] and 0C2 used in the current 
code are 8 and 4. By referring to equations (23)-(25), one 
may see that this means that the SWL limiting process is 
only used in a very few places in the code where 
extremely strong waves arc present. We note here that if 
the extrapolation/interpolations are second order rather 
than third order, the P values of equations (23)— (25) 
would be examined at only two cell boundaries, instead of 
the three as in the above discussion. The extrapolation/ 
interpolation and limiting procedures are essentially the 


same at zone boundaries as at cell boundaries as described 
above. 

3.3 Flux Calculations 

With the left- and right-side values of the primitive vari- 
ables at the cell or zone boundary determined as described 
above, the fluxes across the boundary are determined 
using a Riemann solver, as for any Godunov-type method. 
Figure 5 shows typical Riemann solutions. Zones 1 and 4 
in each case are the left- and right-side conditions, 
extrapolated/interpolatcd as described above. In fig- 
ure 5(a) is shown the usual Riemann solution, with two 
waves, in this case an expansion wave (E) and a shock 
wave (S) and an interface. Two new regions, 2 and 3, are 
created by the Riemann solution dynamics. 

Depending upon the initial conditions in zones 1 and 4, 
the two waves can be expansion and shock (E,S) (as 
shown in fig. 5(a)) or (S,E), (S,S), or (E,E). If two very 
powerful expansions are created in the Riemann problem, 
it may be necessary to include a fifth, minimum pressure, 
zone as shown in figure 5(b). This would be a vacuum 
region for a situation where at least one of the regions 1 
and 4 was a gas or a region at the tensile yield stress if 
both regions l and 4 were dense media. Our code docs not 
have the ability to model spall, that is, to create new free 
surfaces in solids under high-tensile stresses. Rather, the 
solid is modelled as capable of stretching indefinitely at 
the yield stress. While this obviously makes the code 
unsuited for modelling spall situations, both CFD work 
and experimental observations of conditions of the piston 
and projectile in light gas guns strongly suggest that this 
code limitation is of little consequence for almost all light 
gas gun operating conditions modelled. The solutions 
obtained can, of course, he monitored for the occurrence 
of stresses at the tensile yield level in solid material zones. 

The essential components of the Riemann solver are 
expansion and shock solvers. The differential form for a 
general isentropic expansion or compression is given in 
equations (26)-(28) below. 


dp = c^dp (26) 

de = -tdp (27) 

P 

du=±^ (28) 

pc 





Minimum 

pressure 



Figure 5. Riemann problems at cell or zone boundary, (a) Shows the usual problem , with two waves and an interface; 
(b) shows the special problem, with two very strong expansions and a minimum pressure or vacuum region between the 
expansions. 
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For a weak wave (either expansion or compression), with 
Idu/cl not more than 0.03, the above expressions are used 
in the code with little modifications. For stronger waves, 
more accurate expressions are required. 

The calculation of stronger expansion waves starts with 
the determination of the exponent relating the sound speed 
variation with the density variation, i.e., a in equa- 
tion (29) below. 


c 



C| 


KPi) 


(29) 


This is determined by varying the density a small amount 
from the starting condition. (The starting condition is 
denoted by subscript 1 in eq. (29)). It turns out that this 
exponent is, in fact, very nearly constant for a wide range 
of real media even for very deep expansions. Hence, with 
a determined for small density variations and then 
assumed constant for substantial density variations, rather 
good exact integrations of equations (26)-(28) can be 
obtained. These can be readily derived and hence, are not 
given here. These integrations were used for expansions 
in ourRiemann solvers. The integrations can easily be 
carried out to either a specified velocity or pressure. The 
calculated values of the density and energy at the end of 
the expansion are limited, in that they are not allowed to 
go below 1 percent of the corresponding initial values. 
Such limiting need be applied only for a very small num- 
ber of cases in a typical solution. 

Three shock wave solvers have been used with success in 
the code. The first is a linear solver with density limiting; 
equations (30)-(34) are used in this case. 


P2 =Pl + Pl c l| u 2“ u l 


P2 “ Pi 
P2=Pl + — 2 ~ 


PLIM =r pPl 


P2f = min(p 2 ,PLIM) 


e 2=ei+-(pi+P2) 


iPl P2f j 


(30) 

(31) 

(32) 

(33) 

(34) 


The conditions with subscripts 1 and 2 are before and 
after the shock, respectively. Equations (30) and (31) are 
derived directly from equations (26) and (28) above. 

Since the density ratio is known to be limited for strong 
shocks, the derived value of p2 is limited to p2f using 
equations (32) and (33). The values of rp are constants for 
each media input into the CFD code. Reasonable values 
are selected after considerable study of the behavior of 
each media. For gunpowder/gunpowder gas, polyethylene 
and hydrogen, the values of p r chosen were 9.55, 2.2, and 
9.0, respectively. Equation (34) for the energy, e2, is exact 
if the values of p2 and p2f were exact. 

The second solver, a quadratic solver, is presented next. 
This solver was derived independently (ref. 10) by one of 
the authors (DWB), but afterwards was found to also be 
used by workers at Los Alamos (ref. 33). The pressure 



(35) 


P2 =Pl + Pl c l | u 2 _u l| + Ap,(u 2 — u,) 


2 


(36) 


V= I U2 ~ U| I 


(37) 


P2 = 


Pi 


1 i — 

1 + Avj/ 


(38) 


velocity relation of the linear solver, equation (30) is 
modified with the addition of a quadratic term in the 
velocity difference to yield equation (36). With the pres- 
sure and velocity behind the shock known, the continuity 
and momentum equations can then be used to obtain the 
density behind the shock; the result is given in equa- 
tions (37) and (38). The energy behind the shock is 
obtained from equation (34), as was done for the linear 
solver. 

The third solver option is our “exact” shock solver, which 
iterates to obtain the post-shock conditions. Using the 
“exact” solver option, when \\f is less than 0.03, the code 
defaults to the linear solver. The “exact” solver guesses 
p2, limited by the values p \ and rpp j and uses the conti- 
nuity, momentum, and energy equations to calculate P2, 
U2, and e2- The EOS is then used to generate p'2(P2> e 2)- 
P2 is then iterated until p2 equals p’2. The iteration is 
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stopped when the number of iterations reaches 12 or the 
normalized P2 error reaches 0.005, whichever occurs first. 
Essentially no improvement was found during code 
development by increasing the number of iterations to 24 
and decreasing the acceptable p2 error to 0.001 . Thus, the 
accuracy given here is believed to be quite sufficient for a 
shock wave solver used within a Riemann solver. 

A set of code development checkout runs were made with 
the impact of two polyethylene slabs at a closing velocity 
of 20 km/sec. One of the slabs is backed by an infinitely 
stiff plate. The outside surface of the other slab is free. 
Thus, the problem produces two initial shock waves, one 
reflected shock wave and a reflected expansion wave. 
Exact solutions can readily be calculated for this problem 
up to a certain point in time. A large improvement in the 
solutions was noticed when shifting from the linear to the 
quadratic shock solvers described above. A smaller, but 
still noticeable improvement in the solution was noticed 
when shifting from the quadratic to the “exact” shock 
solver. Hence, in the solutions presented herein, the 
“exact” shock solver was always used. There is a penalty 
in CPU time in shifting towards the better shock solvers. 
For the slab impact code checkout runs, the run times 
were found to be essentially in the ratio 1 :2:4 for use of 
the linear, quadratic, and “exact” shock solvers, 
respectively. 

We now move on to the method used to obtain a complete 
Riemann solution, with all 4 (or 5) zones shown in fig- 
ure 5. First, complete expansion waves are calculated 
from the two initial zones 1 and 4 to see if the situation 
with the minimum pressure zone of figure 5(b) exists. If it 
does, then, at that point, the complete solution with all 
waves and regions has already been found. In the vast 
majority of cases, the minimum pressure zone is not 
required and we have the situation shown in figure 5(a) 
with four zones (though not, of course, necessarily with 
the particular pair of waves — expansion/shock — shown in 
the figure). The following procedure is then used to solve 
the problem. Two linear (acoustic) p-u equations (similar 
to eq. (30)) are written across the two waves between 
regions 1 and 2 and also between regions 3 and 4. These 
two equations are solved by setting U2 = U3 and p2 = P3 to 
yield a value for U2 - U3. (The p2 = P3 values are 
discarded.) Then, the higher-order expansion and shock 
solvers discussed above are used from u\ to U2 and also 
from U4 to U3 to generate new values for p2 and P3, which 
do not, in general, agree with each other initially. The key 
to our procedure is the technique of selecting the updated 
value of U2 = U3 to be used in the next iteration. This is 
done by again solving a Riemann problem using the 
acoustic p-u relations, but starting now from the previous 
P2,U2 and P34J3 conditions found using the higher-order 
shock and expansion solvers. This leads to a new value 


for U2 — U3. Once this new value of U2 = U3 is found, the 
previous solutions for states 2 and 3 are completely 
discarded and a new higher order set of calculations of 
states 2 and 3 is made to the updated U2 = U3 values. This 
procedure has been found to be exceedingly robust and to 
converge very rapidly. The procedure is repeated until the 
normalized error between p2 and p3 is less than 0.001 or 
the iteration number is eight, whichever comes first. In 
most cases, the procedure was found to converge in two 
or three iterations. Again, increasing the iteration number 
and improving the accuracy further were found to make 
no perceptible difference to the solutions. 

With the complete Riemann solution in hand, it now 
remains to select the region in which the cell boundary 
lies to calculate the flux across said boundary. Since the 
cell and zone boundaries slide, in general, the line fol- 
lowed by the boundary in question will usually not be ver- 
tical in figure 3. (We will refer to this line as the “world 
line” of the boundary in question.) If the world line of the 
boundary is within one of the zones l, 2, 3, 4, or 5, the 
flux is simply calculated using the primitive variables for 
that zones. Component mass fractions (and transverse 
velocities, if a multidimensional solution were being car- 
ried out) arc simply those of region 1 or 4, depending on 
whether the boundary world line is to the left or right of 
the interface. 

If the boundary world line lies within an expansion wave 
system, an additional expansion integration is carried out 
exactly to the world line to generate a new set of primitive 
variables which arc then used to calculate the fluxes. The 
other option would be not to perform the extra integration 
and to assign the boundary world line to one or the other 
of the zones bounding the expansion wave depending on 
whether the world line lay to the left or right of some 
defined “middle” of the expansion wave. The fluxes 
would then be calculated directly from the primitive vari- 
ables in the chosen zone. We have found virtually no dif- 
ference in the code performance between these two 
techniques. Nevertheless, we have continued to use the 
technique with the extra integration, since it is, in princi- 
ple, more correct and adds only a fraction of one percent 
to the CPU time requirements of the code. 

3.4 Boundary Conditions 

Calculation of the boundary fluxes comprises two compo- 
nents. First, at each end of each zone, conditions (i.e., the 
state vector, U) are calculated in a “ghost cell” which lies 
just outside the zone. A brief discussion of the use of 
so-called “ghost” cells is appropriate at this point. At both 
ends of a zone, just outside the zone, one or more ghost 
cells are added to the chain of real cells which form the 
interior of the zone. The cell-center values of the ghost 
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cells arc chosen so that a Riemann solution between the 
real and ghost cells straddling the zone-boundary repro- 
duces the proper zone-boundary condition. The use of 
ghost cells is extremely convenient when higher-order 
interpolalions/extrapolations must be made near the zone 
boundary, since these can now be done essentially in the 
same way as is done in the interior regions of the zones. 
We now return to the main-boundary condition 
discussion. 

Using the U values in the ghost cells, as well as those in 
the real cells of the zonc(s), extrapolations/interpolations 
arc made to the zone boundary. The fluxes there are calcu- 
lated using Riemann solver techniques analogous to those 
used for the cell boundaries within a zone, as described in 
section 3.3. There is a certain degree of freedom in choos- 
ing the ghost cell values in a number of cases and this 
freedom becomes much greater when two ghost cells are 
used. For this reason, we have limited the code to the use 
of only one ghost cell and hence, second-order 
cxtrapolalions/interpolations to the zone boundaries are 
used under some conditions. In addition, for some be used 
only for code check-out, we have used cell-center (i.c., 
first order) values to calculate the zone-boundary fluxes. 
Below, we will briefly go through the be used in the code. 

For the infinitely stiff wall at zero velocity be (be #1 ), the 
ghost cell U (slate variable) is created by taking the U of 
the last, adjacent real cell in the zone and simply reversing 
the velocity, i. c., replacing u by -u. For the point mass 
projectile be (be #10), the ghost cell U is created by taking 
the U of the last, adjacent real cell in the zone and replac- 
ing u by 2u pm j-u, where Up ro j is the projectile velocity. 
For the free surface be (be #2), the ghost cell U is created 
by solving a Riemann problem expanding to zero pres- 
sure. Solutions of all Riemann problems involved in 
determining the be are done using the same techniques 
already described in section 3.3. There are three bes are 
used for internal zone boundaries — moving boundary, 
different media (be # 4); fixed boundary, same media (be 
#5); and moving boundary, same media (be #6). In all of 
these cases, the ghost cell Us arc created by solving a 
Reimann problem between the two adjacent real end cells 
of the two zones. For some of the preceding bes, some of 
the extrapolations/interpolations subsequently used to the 
boundary in question to calculate the boundary fluxes arc 
of second order, rather than third order. 

Four other bes were used for code development and 
debugging work. The nonreflccting outflow condition (be 
#3) sets the ghost cell U equal to that of the adjacent real 
cell. The supersonic inflow be (be #7) keeps the ghost cell 
U at a fixed value specified at the beginning of the CFD 
solution. In the subsonic inflow be (be #8), the ghost cell 
pu and H = c + p/p + u 2 /2 are kept at specified, fixed val- 


ues and a Riemann solution is performed between the 
ghost cell and the adjacent real cell to determine com- 
pletely the ghost cell U. Similarly, in the subsonic outflow 
be (be #9), the ghost cell p is kept at a specified, fixed 
value and a Riemann solution is performed between the 
ghost cell and the adjacent real cell to determine com- 
pletely the ghost cell U. For all of the preceding four bes, 
the extrapolations/interpolations subsequently used to the 
boundary in question to calculate the boundary fluxes arc 
first-order (i.c., cell center) values. 

For the blind end of the breech, be #1 is used. For both 
sides of the piston, be #4 is used. For the main 
diaphragm/projectile base condition, be #1 is used before 
rupture. After rupture, be #10 is used for point mass pro- 
jectiles and be #4 for projectiles with an internal zone and 
celling. For the free surface of an internally zoned projec- 
tile, be #2 is used. If there is an extra diaphragm in the 
pump tube, be #1 is used there before break and be #6 
after break. 

3.5 Multiple Zoning 

The computational domain is divided into several zones. 
Regions of different media, i.e., gunpowder (plus powder 
gas), polyethylene piston material, hydrogen pump-tube 
gas are always kept in different zones, to avoid the EOS 
difficulties that would occur if two different media 
co-existed in one cell. This problem is avoided by sliding 
the grids so that the grid zone boundaries move with the 
media boundaries. In some cases, more than one zone can 
be used for a single media, in particular, with the pump- 
tube gas. The gridding zones must, of course, change size 
in “accordion’’ fashion to follow the expansions and con- 
tractions of the media zones as the gun cycle progesscs. 

The cell sizes are arranged so that across a zone boundary, 
the zone end cell sizes are equal. This is done as follows. 
Every second zone has its cell sizes calculated at each 
time step independently of all other zones. Then the 
remaining zones have their end cell sizes calculated to 
match the end cell sizes (already calculated) immediately 
across the zone boundaries. Three methods of calculating 
the cell sizes are available for the first set of zones 
(calculated independently of all other zones). The first 
two are: (1) uniform cell sizes and (2) cells increasing (or 
decreasing) in size by a given, constant ratio from one end 
to the other of the zone. Method (3) is a combination 
technique. A critical cell size is chosen. If the zone length 
divided by the number of cells is less than this size, uni- 
form cell sizing is chosen. If it is greater than this length, 
the cell size at a chosen end of this zone is set to the criti- 
cal size and, moving away from this cell, the cell sizes 
successively increase by a ratio chosen to exactly fill the 
zone with cells. This ratio varies with every timestep. The 
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purpose of this complex zoning technique is (a) to allow 
such a zone to keep a reasonably small cell size (to obtain 
adequate resolution) at one end when the zone is greatly 
expanded and (b) to avoid extremely tiny cells and the 
consequent great reduction in code timestep when the 
zone is greatly compressed. 

The remaining zones are calculated either by method 
(4) or method (5), as follows. In method (4), if the zone is 
an end zone, it is calculated with geometrically increasing 
(or decreasing) cell sizes from one end to the other (con- 
stant size ratio between adjacent cells). By adjusting the 
size ratio between adjacent cells, the cell size at the end of 
the zone abutting the neighbor zone can be adjusted to 
match that of the end cell of that zone. This size ratio 
varies with every timestep. If the zone is between two 
other zones, the previous technique does not provide 
sufficient degrees of freedom to match the end cell sizes 
across both zone boundaries. This problem is resolved by 
using method (5) — zoning the region with two different 
subregions each with geometrically varying cell sizes with 
two different size ratios between adjacent cells for the two 
regions. With two different size ratios available to vary, 
one can match the end ceil sizes across each zone 
boundary of the zone and across the boundary between 
the two subregions. 

The type of zoning must be selected by the program 
operator. For a three-zone (gunpowder, piston, and hydro- 
gen) two-stage gun solution, good results have been found 
with the following zoning methods: gunpowder, method 
(4); piston, method ( i ); and hydrogen, method (4). For a 
four-zone (gunpowder, piston, hydrogen #1, and hydro- 
gen #2) two-stage gun solution with an extra diaphragm in 
the pump tube between the two hydrogen zones (ref. 1 1) 
good results have been found with the following zoning 
methods: gunpowder, method (4); piston, method (1); 
hydrogen #1, method (5); and hydrogen #2, method (3). 

3.6 Piston Jam Versus Bounce 

The amount of friction calculated in the CFD code on the 
piston in generally not enough to prevent the piston from 
bouncing back after maximum gas compression. For our 
operating conditions of the NASA Ames 1.5 in. gun, it is 
observed that the piston, in fact, jams in the high-pressure 
coupling rather than bounces. To account for this, an 
option is available in the code which changes the be at the 
piston front to the infinitely stiff wall at zero velocity be 
#1 in section 2.7 (on both sides of the zone boundary) 
when the velocity of the piston front reaches zero. Most of 
the two-stage gun solutions presented herein use this 
option. It will be shown that this option produces excel- 
lent agreement between theory and experiment. The 
behavior of the piston exposed to high pressures appears. 


as will be discussed at a later point, to be very poorly 
understood and to disagree very substantially with low 
speed measurements of the yield stress and Young’s 
modulus of the piston material. 

4 Code Validation Versus Analytical 
Solutions 

The code validation consists of two parts — validation 
against analytical solutions and validation against data 
from actual two-stage gun firings. In this section, we will 
outline the validations performed against given analytical 
solutions. 

4.1 Riemann’s Shock-Tube Problem 

The first-validation work is against analytical solutions 
for Riemann’s shock-tube problem. Figure 6 shows (a) the 
initial configuration with the high-pressure driver tube, 
the low-pressure driven tube and the diaphragm and (b), 
the x-t wave diaphragm. Upon rupture of the diaphragm, 
an expansion wave system moves into the driver-tube gas 
and the incident shock wave, followed by the driver- 
driven gas interface moves into the driven tube. The inci- 
dent shock wave subsequently reflects from the driven- 
tube end wall, as shown. As far as is shown in figure 6, 
the wave processes can be modelled exactly analytically 
and thus serve as an excellent test for the CFD code. 

The Riemann problem modelled analytically and by the 
CFD code used ideal gases, included no wall-friction or 
heat-transfer effects and had the following parameters: 

Molecular weight = 2 

Specific heat ratio = 1 .4 

Temperature = 700 K 

Pressure = 6.898 x 10 8 dynes/cm 2 

Length =2100 cm 

Driven tube/driven gas - 

Molecular weight = 29 

Specific heat ratio = 1 .4 

Temperature = 295 K 

Pressure = 8.669 x 105 dynes/cm 2 (or 

fraction thereof) 

Length = 2600 cm 

With Riemann’s shock tube problem, we can compare the 
analytical and CFD code solutions for the (1) incident 
shock, (2) the driver rarefaction, and (3) the reflected 
shock. 
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Figure 6. Riem arm's shock-tube problem, (a) Shows initial configuration with driver and driven tubes and diaphragm ; 
(b) shows x-t wave diagram. 


CFD solutions were carried out with 20 cells each in the 
driver and driven tubes for the above pressure ratios (Rp) 
of (approximately) 10^, and also, by reducing the driven 
tube pressure by successive factors of 10, for pressure 
ratios of (approximately) I0 4 , 10^, and 10*\ The solutions 
were then compared with exact analytical results. The 
CFD solutions were excellent at Rp = 10^, showed modest 
deterioration at R p = 10 4 and progressively more serious 
deterioration at R p = 10 5 and 10 6 . As an example of dete- 
rioration, the peak-to-peak amplitude of the pressure 
oscillations behind the incident shock in the Riemann 
shock tube problem was 1-2 percent for R p = l(P 
and 10 4 , 10 percent for R p = 10^ and 40 percent for 
Rp = 10 6 . Nevertheless, the code was robust in that it pro- 
duced solutions at R p up to 10^ with only 20 cells per 
zone, and that the degradation of the solutions as R p 
increased from 10 4 to 10 6 took place in a smooth, consis- 
tent fashion. The number of cells per zone was then 
increased to 40 and to 80. With 80 cells per zone, the 
solutions were now excellent for R p of 10 4 , and the rather 
large entropy (temperature) errors of the solutions for R p 
of 10 5 and 10 6 with 20 cells per zone were much reduced. 
For example, (spurious) oscillations in temperature in the 
reflected shock region for R p = 10 6 were reduced from 


67 percent to 17 percent when the number of cells was 
increased from 20 to 80. At the same time, the oscillations 
behind the incident shock referred to earlier were reduced 
from 40 percent to 12 percent (for R p = 10 6 ). Finally, we 
switched from uniform cells in the driven gas to having 
these cells decreasing in size by a constant ratio from the 
end wall to the interface. With an overall decrease in cell 
size within the zone by a factor of 36, still further 
improvements in the solution for R p = 10^ were obtained. 
The CFD code issues discussed above regarding the 
Riemann shock tube problem are very similar to those 
presented in reference 34 for solutions with between 200 
and 3200 cells (total). 

4.2 Plate-Slap Problem 

In this problem we consider two plates of polyethylene 
impacting each other at a closing velocity of 20 km/sec. 
One plate (the anvil) is initially at rest and is backed by an 
infinitely stiff wall. The rear surface of the moving plate 
(flyer) is completely free. Figure 7 shows (a) the initial 
configuration with the flyer, anvil, and infinitely stiff wall 
(b), the x-l wave diaphragm. Upon impact of the flyer and 
the anvil, shock waves are created in both. The anvil 
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Figure 7. Plate-slap problem, (a) Shows initial configuration with flyer, anvil, and infinitely stiff wal;, (b) shows x-t wave 
diagram. 


shock wave reflects as a shock wave from the infinitely 
stiff wall and the flyer shock wave reflects as a complete 
centered rarefaction wave system from the flyer free sur- 
face. Again, as far as is shown in the figure, the w^ave pro- 
cesses can be modelled exactly analytically and thus serve 
as an excellent test for the CFD code. 

The dense media EOS #3 (see sec. 2.2) is used for the 
polyethylene. This test case allows study of (I) the initial 
pair of shock waves in both plates, (2) a shock wave 
reflected from the infinitely stiff wall, and (3) a complete 


expansion wave reflected from the free surface. Compar- 
isons were again made between the CFD code solutions 
and exact analytical solutions. Solutions were carried out 
with a total of 40, 80, and 160 cells. The solutions were 
excellent, especially for the cases with 80 and 160 cells. 
Some small wiggles, of the order of 2-3 percent, were 
observed immediately downstream of shock waves, but 
quickly damped out, and some entropy errors 
(5-10 percent maximum density error) were created 
where the shock waves were created or reflected. 
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Otherwise, the solutions were very accurate, with all 
shock- and expansion-wave velocities and pressures very 
well predicted. 

4.3 Convergent-Divergent Nozzle 

In the preceding code validations, the cell area has been 
constant, i.e., the solution has been completely one- 
dimensional. To test the code with area changes, we 
examine steady, frictionless flow of ideal air through a 
convergent-divergent nozzle. The nozzle area decreases 
by a factor of 16 to the throat and then increases by a fac- 
tor of 16 back to its original area. The total length of the 
portion of the nozzle with area variations is 32 cm and the 
nozzle diameter varies from 10 cm to 2.5 cm and then 
back to 10 cm. The CFD code calculation was done with 
90 cells filled with air. Both air and cells are pushed 
through the nozzle by pistons located upstream and down- 
stream of the nozzle. The air does not move with the same 
velocity as the cells, since this is not a Lagrangian code. 
The initial cell conditions are set to the exact quasi-one- 
dimensional solution and then about 40 cells are pushed 
through the nozzle. The solution is re-examined against 
the exact quasi-one-dimensional solution to assess code 
accuracy. (Some of the cells make a complete traverse of 
the portion of the nozzle with area changes.) The ideal 
pressure ratio across the nozzle is about 300 to 1 and the 
exit Mach number is about 4.5. 

The errors were assessed primarily by evaluating the 
entropy and enthalpy throughout the nozzle. The maxi- 
mum enthalpy (or temperature) error was 0.5 percent and 
the maximum entropy error was AS/R = 0.02, which 
would correspond to a pressure error of 2 percent. This 
was judged very satisfactory for a nozzle with a pressure 
ratio of 300. With half as many cells, the errors are con- 
siderably larger, with a maximum enthalpy error of 
2.7 percent and a maximum entropy error of AS/R = 0. 1 7. 

4.4 Gunpowder Burn 

Gunpowder burn in the CFD code within a closed bomb 
with 10 cells was compared with results obtained from a 
very highly accurate (small timestep) direct (non-CFD) 
integration of the powder-burn equations. In both calcula- 
tions, the shape of the powder grain changes as it is con- 
sumed. All results agreed within 0.5 percent between the 
two sets of calculations. 

4.5 Friction and Heat-Transfer Calculations 

No simple analytical solutions are available against which 
to test the friction and heat-transfer calculations. How- 


ever, special checks were made of the friction and heat- 
transfer calculations for gases and dense media and of the 
nonequilibrium turbulent kinetic energy calculation for 
gases, as follows. For each of these calculation tech- 
niques, for one cell, for one timestep, a complete hand 
calculation was made to six decimal places and the results 
compared with extensive special diagnostic print-outs 
from the CFD code. A few errors were found in this way 
and were corrected. 

5 Code Validation Versus Actual Gun Data 

Sections 5. 1-5.5 discuss mainly code validation with 
experimental data from actual two-stage gun firings. Sec- 
tion 5.6 discusses mainly important predictions made by 
the code which we cannot, at present verify experimen- 
tally, such as the pressure histories in the high-pressure 
section of the gun and at the projectile base. Finally, sec- 
tion 5.7 discusses apparent anomalous yielding and 
deflection behavior of the pump-tube piston. 

5.1 Gun Configuration and Operating Conditions 

Figure 8 is a sketch (not to scale) showing the dimensions 
of the Ames 1.5 in. gun (with the benchmark piston) as 
modelled in the CFD code. The high-pressure coupling 
contraction is actually modelled with two cones with a 
slight break in the slope at the 3358.58 cm station at a 
diameter of 1 1 .35 cm. This is not shown in figure 8 for 
clarity. However, using two cones models the actual mea- 
sured contraction dimensions considerably better than if 
only one cone is used. Table 1 gives the benchmark oper- 
ating conditions of the gun. 

Figure 9 shows (a) the real piston and (b) the one- 
dimensional piston used in the CFD code. 

The one-dimensional piston has the same land lengths and 
land and shank diameters as the real piston and its length 
has been adjusted so that its mass is equal to that of the 
real piston. 

Table 2 summarizes the experimental and CFD data for a 
number of shots of the Ames 1.5 in. gun. The rupture 
pressures of the break valves (diaphragms) were calcu- 
lated according to reference 35. Shots 600-607 (omitting 
shot 604) were at nominally identical gun operating 
conditions. The data from these shots were averaged to 
produce an “average” shot listed as “AV. 600-607” in the 
table. The powder-energy release, the powder-bum rate 
and the piston friction were adjusted to match the 
observed powder-chamber pressure histories and the pis- 
ton velocity. 
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Figure 8. Gun modelled in CFD study. Not to scale. All dimensions are in centimeters. Station numbers are distances 
from blind end of breech. DIA denotes diaphragm (break valve). 


Table I. Ames i.5 in. benchmark gun operating conditions 


Powder 

Hercules HC-33FS 

Powder mass 

3000 gm 

Piston material 

High-density polyethylene 

Piston mass 

21,420 gm 

Pump-tube hydrogen pressure 

3.104 bar 

Break- valve break pressure 

1216 bar 

Sabot material 

Lexan 

Total launch mass 

29.7 gm 



Figure 9. Representative pump tube piston for NASA Ames 1.5 in. gun (not to scale). Dimensions in centimeters, 
(a) Actual piston; (b) one-dimensional piston for CFD code. 
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Tabic 2. Experimental and CFD data for shots on NASA Ames 1 .5 in. gun 


Shot no. 

Powder 

mass 

(gm) 

Piston 

mass 

(gm) 

Project. 

mass 

(gm) 

Break 

valve 

press. 

(bar) 

Hydrogen 

pressure 

(bar) 

CFD 

piston 

velocity 

(m/s) 

Exper 

piston 

velocity 

(m/s) 

CFD 

projectile 

velocity 

(km/s) 

Exper. 

project. 

velocity 

(km/s) 

600 

3000 


30 

1216 

3.104 





601 

3000 

21,420 

30.13 

1216 

3.104 


683.2 


6.765 

602 

3000 

21,420 

30.14 

1216 

3.104 


685.3 



603 

3000 

21,420 

30.3 

1216 

3.104 


685.3 


6.887 

605 

3000 

21,420 

29.81 

1216 

3.104 


689.5 


6.638 

606 

3000 

21,420 

29.68 

1216 

3.104 


685.9 



607 

3000 

21,420 

29.72 

1216 

3.104 


685 


6.714 

AV. 600-607 

3000 

21,420 


1216 

3.104 

685.8 

685.6 

6.744 

6.722 

613 

2835 

17,064 

29.015 

1216 

3.104 

742.8 

735.9 

7.002 

6.818 

614 

2835 

17,064 

29.317 

1216 

3.104 

742.8 

735 

7.002 

6.829 

615 

2810 

17,064 

29.987 

690 

3.104 

739.4 

728.1 

6.951 

6.307 

616 

2835 

17,064 

29.57 

813 

4.711 

710.3 

716.1 

7.092 

7.217 

618 

2835 

17,064 

33.579 

785 

4.7! I 

710.3 

721.8 

7.064 

7.059 


Table 3. Selection of powder and piston parameters 


Gun 

Ames 1 .5 in. 

Ames 0.28 in. 

Powder 



Maker 

Hercules 

Du Pont/IMR 

Type 

HC-33FS 

4227 

Energy factor 

0.97 

0.97 

Burn rate factor 

1.891 

1.581 

Piston friction model 



Lands included 

No 

Yes 

Land diameter 

N/A 

0.235 x actual + 0.765 x bore 

Friction coeff. factor 

0.396 

1 


5.2 Selection of Powder and Piston Parameters 

Table 3 shows the parameters chosen for the Ames 1.5 in. 
and 0.28 in. guns. In the “powder” section of the table, the 
first-two lines give the maker and the type of the powder. 
The next line gives the energy factor — this is the factor by 
which the powder energy release calculated from the 
manufacturer’s data (see sec. 2.6) must be multiplied to fit 
the experimental data. Next is given the burn rate factor — 
this is the factor by which the powder burn rate given by 
the manufacturer (see sec. 2.6) must be multiplied to fit 
the experimental data. 

As mentioned above, the powder energy release, the 
powder burn rate and the piston friction were adjusted to 
match the observed powder chamber pressure histories 
and the piston velocity. Matching of the powder-chamber 


pressure history was made to the pressure rise and to the 
pressure fall from maximum to about 1/3 of the maximum 
value. The adjustment of these parameters requires the 
insertion of “fixes” into the code and we will now outline 
the reasons why we believe these are necessary. The 
adjustment of the powder energy release is very small 
(3 percent) and is the same for both types of powder and 
may possibly represent incomplete burning of the powder. 
It is also necessary to point out that fits to the experimen- 
tal data that are almost as good as those obtained here are 
probably obtainable with no adjustment to the powder 
energy release. It is suggested that the powder energy 
factor be taken to be between 0.97 and 1 .0. 

On the other hand, the burn rate parameter is much 
greater than unity (1 .58-1 .89) and is quite different for the 
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two powders studied. There is a feedback mechanism in 
that higher burn rates raise the powder-chamber pressure 
which, in turn, causes still further increases in the bum 
rate. Hence, a burn rate parameter increase from 1 .0 to, 
say, 1.7 produces a much larger relative increase in cham- 
ber pressure. The upshot of this is, that if one uses the 
burn rates given by the powder manufacturers in the code, 
the chamber pressure bears almost no resemblance to that 
actually observed, yielding a peak pressure only half (or 
less) than that observed and a pressure history very unre- 
alistically stretched out in time. 

Two effects likely contribute to the increased burn rate 
(over that predicted using the powder manufacturers’ 
specified constants in equation (10) in sec. 2.6) in our 
two-stage gun breeches. First, the maximum pressures in 
the breeches of Ames’ two-stage guns range from 
400-1000 bar, whereas in typical military gun systems, 
the maximum pressures are often 3000 bar or higher 
(ref. 36). The manufacturers’ bum-rate parameters are 
obtained from closed bomb or strand-burner tests and they 
are usually made (ref. 30) to emphasize the pressure range 
of 600 to 2000-2700 bar which is more applicable to mili- 
tary weapons. Some of the gun propellants for which data 
is given in references 30 and 37 show anomalies in the 
pressure range of 400-1000 bar which would translate to 
higher burn rates at these pressures than those which 
would be predicted by the powder manufacturers’ equa- 
tions. This point is also discussed in reference 38. Stated 
another way, it is possible that the manufacturers’ burn- 
rate expressions may not be very applicable to the rela- 
tively low breech pressures of our two-stage guns. 

A second effect is that in guns, there is gas flow over the 
propellant grains and the resulting increases in convective 
heat transfer cause corresponding increases in the burn 
rate (“erosive burning”) over those which would be mea- 
sured in closed vessels and are the usual basis of the 
manufacturers’ bum rate equations. This point is dis- 
cussed further in reference 39. The benchmark bum rate 
in our code is the standard expression given by the pow- 
der manufacturer. Further, the code, by choice, does not 
include the calculation of any erosive burning effects. The 
burn rate must be increased substantially above that given 
by the powder manufacturer in order to match the actual 
observed powder-chamber pressure histories. This is 
likely a result of the two effects discussed above. The fac- 
tors of rate increase (“burn rate factor”) observed to pro- 
duce good agreement between theory and experiment for 
the two powders used in the Ames 1.5 in. and 0.28 in. 
guns (see table 3) were 1.89 and 1.58. Based on these 
numbers, it is recommended, in fitting a new powder in a 
new gun using our code, that one start with a burn rate 
factor of about 1 .75 and iterate in to make the experimen- 


tal and CFD powder-chamber pressure histories and the 
piston velocities agree. 

In addition, the piston motion in the pump tube was 
found, in general, to take place as though the piston fric- 
tion was substantially less than that calculated as dis- 
cussed in section 2.5 and Appendix C. The real piston 
behavior can be much more accurately modelled by 
adjusting parameters in the piston model. One may 

(1) reduce the land diameter to reduce the jamming forces, 

(2) remove the lands entirely, or (3) reduce the coefficient 
of friction of the piston material against the steel pump 
tube. As indicated in table 3, for the 1.5 in. gun piston, the 
lands were removed entirely and the coefficient of friction 
was reduced to 0.396 times the normal values which 
would be used according to the models of section 2.7 and 
Appendix C. For the 0.28 in. gun, the land diameter was 
reduced as indicated in the table, but the friction coeffi- 
cient was not changed from its normal value. [One should 
note here that the (room temperature) land diameter is 
greater than the tube bore and the pistons are inserted into 
the pump tube after having been cooled in a freezer and 
then are allowed to come up to room temperature and 
expand and jam in the pump tube.] The two methods cho- 
sen to reduce the piston friction shown in table 3 differ for 
historical reasons only — there is no reason to prefer one 
over the other. 

There are a number of major difficulties in modelling the 
piston friction. First, the piston land surface may melt, 
reducing the friction coefficient below that used in the 
model (sec. 2.5 and App. C). Second, the lands will wear 
away progressively, thereby reducing the land jam in the 
pump tube, over the course of the piston stroke. Both of 
these phenomena are very difficult to model. Third, there 
is considerable evidence (see sec. 5.7) that the high strain 
rate deflection and yielding behavior of the piston mate- 
rial is very different from that predicted using material 
parameters measured in low strain rate situations. These 
latter parameters are used in our model (following ref. 4), 
lacking any data on high strain rate deflection and yield- 
ing behavior of the piston material. In general, the exper- 
imentally observed behavior of the piston indicates a sub- 
stantially lower friction drag force than that which would 
be calculated applying the model of section 2.5 and 
Appendix C, directly, without modification. For the two 
guns modelled, the piston of the larger gun requires a 
much greater modification (reduction) in the calculated 
piston friction force to match the experimentally observed 
data. This may be due to the greater relative precision and 
smoothness of both the piston surfaces and the pump tube 
bore. This would tend to occur since the surface finishes 
and manufacturing tolerances tend to be nearly constant, 
regardless of component size. Hence, for a gun bore 
which is about five times larger (the Ames 1.5 in. gun), 
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the surface finishes and tolerances, normalized to the gun 
dimensions, are about five times better. 

To reduce the piston friction to obtain agreement between 
CFD and experimental data for a new gun using our code, 
it is suggested to reduce the piston land diameter and the 
piston material friction coefficient as discussed above and 
shown in table 3. Since there are wide differences in the 
amounts of friction reduction required to obtain a good 
match between CFD and experiment and since piston 
friction is so poorly understood, we do not recommend 
even a first guess of how much to reduce the friction coef- 
ficient to model a new gun. This must be determined on a 
case by case basis. We point out that in determining the 
burn-rate factor and the piston friction reduction to be 
used with our code for a new gun, both parameters must 
be varied to get good agreement between the CFD and 
experimental chamber pressure histories and piston veloc- 
ities. Typically, about 10 runs must be made with well 
chosen parameters to obtain good matches. It is suggested 
that a single, well documented, gun operating condition be 
chosen against which to adjust the burn rate factor and the 
piston friction reduction and that thereafter, these parame- 
ters be kept constant (for a given gun and powder type). 

The performance of the gun is very sensitive to the piston 
velocity and hence to the piston friction. Thus, the deter- 
mination of the burn rate factor and the piston friction 
reduction, as discussed above, are very important to 
enable one to obtain good code predictions. It turns out 
that the performance of the gun is very insensitive to the 
method of calculation of projectile friction (within rea- 
sonable limits). We have made CFD runs with nominal 
projectile friction, twice nominal projectile friction and 
zero projectile friction and the muzzle velocities were 
essentially identical for all three cases. This is believed to 
be due in part to the fact that average pressures acting on 
the projectile are much larger relative to the friction forces 
than those acting on the piston. (For much of the piston 
stroke, the pressures on both sides of the piston are quite 
low.) A second effect is as follows. When the projectile 
friction is increased arbitrarily (say, doubled), the first 
part of its trajectory is traversed at a lower velocity, as 
expected. However, it then stays closer to the high- 
pressure gas reservoir for a longer time, permitting the 
pressure waves from the reservoir to catch up to it more 
effectively than if it moved off more rapidly (with lower 
friction). The result is that, with the increased friction, the 
projectile is accelerated more effectively in the later part 
of its in-barrel trajectory. The overall result appears to be 
a near independence of muzzle velocity over a fairly wide 
range of assigned projectile friction values (within 
reasonable limits), as mentioned above. The calculations 
presented herein for the Ames 1 .5 in. gun were made with 


the projectile diameter set exactly equal to the barrel 
diameter, and this is suggested as a starting point when 
applying our code to a new gun configuration. 

5.3. Code Validation with Data from the Ames 1.5 in. 
Gun 

As mentioned in the preceding section, the powder energy 
factor, the powder burn rale factor and the piston friction 
were tuned to make the powder-chamber pressures histo- 
ries and piston velocities, calculated using the CFD code, 
match those for the benchmark gun operating condition, 
which is that of the average of shots 600-607 (see 
table 2). From this table, we see that the CFD piston 
velocity, in fact, is within 0.3 percent of the mean experi- 
mentally observed value, as would be expected, since it 
was tuned to agree. Figures 1 0 and 1 1 show the experi- 
mental and CFD powder-chamber pressure histories. The 
experimental data is for shot 607. Figure 8 shows the pres- 
sure histories out to 22 msec and figure 1 1 shows the his- 
tories out to 62 msec with a logarithmic pressure scale. 

Bit noise becomes progressively more evident in the 
experimental data as the pressure decreases. The 
pressure difference corresponding to one bit is about 
1 x 10 7 dy/cm 2 . The horizontal “lines” in the exper- 
imental data of figure 1 1 are due to switching between 
several bit levels; with the very large number of experi- 
mental data points (15,000), the data points merge into 
apparently solid lines. The matching of the pressure 
histories by the adjustment of the powder parameters and 
piston friction was done only out to a lime of about 
15 msec. The excellent agreement of the pressure histories 
out to pressures an order of magnitude lower and times 
out to 62 msec represent part of the validation of the 
accuracy of the code. 

An important validation of the code is the comparison of 
the experimental and CFD projectile muzzle velocities. 
These are given for the benchmark gun operating condi- 
tion (the line marked “AV. 600-607”) in table 2. The CFD 
and experimental results agree within 0.5 percent. It is 
important to realize that the CFD projectile muzzle veloc- 
ity was not tuned in any way, as was the CFD piston 
velocity. It simply results directly from the CFD analysis. 

The code was also validated by comparing CFD predic- 
tions and experimental measurements in the pump tube 
made at station 3179.1 cm (see fig. 8) in the tube. The 
measurements were made using a PCB model 1 19M39 
pressure transducer and a Kistler model 566 charge ampli- 
fier, as were the measurements of the powder-chamber 
pressure. Figures 12 and 13 show the experimental and 
CFD pump-tube pressure histories. The experimental data 
is for shot 607. Figure 12 shows the pressure histories 
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Figure 10. Experimental and CFD powder-chamber pressure histories. 
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Figure 11. Experimental and CFD powder-chamber pressure histories. 


from 20 to 65 msec with a logarithmic pressure scale and 
figure 13 shows the histories from 48 to 57 msec. Bit 
noise is particularly evident in the lower pressure ranges 
of figure 12. Six shock waves appear in the experimental 
pressure trace up to the time that the piston covers the 
transducer and all six waves are very well predicted by the 
code, with timing errors of perhaps 1-1.5 msec for the 
ilrst wave, 0.5 msec for the next two waves and 


0.2-0. 3 msec for the last three waves. Many details of the 
pressure history are very well captured, such as gradual 
pressure rises between waves 3 and 4 and 4 and 5 
(between 49 and 54.5 msec) and the curved “hooked” 
shape of wave 5 (at about 54.8 msec). 

The experimental pressure levels are consistently about 
15 percent higher than the predictions. It is believed that 
this is a transducer calibration problem, rather than any 
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Figure 12. Experimental and CFD pump-tube pressure histories. 
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Figure 13. Experimental and CFD pump-tube pressure histories. 


code difficulty. For example, if the code were that far off, turn-down ratio, there can be some changes in the trans- 
the timing of the waves would be very much in error and ducer calibration factor and (2) we arc using an uncon- 

other predictions of the code, such as projectile muzzle ventional transducer measuring system, with a PCB 

velocity, would also be substantially in error. The Iran- transducer and a very old Ki slier charge amplifier. It is 

ducer calibration factor could be in error because (1) we possible that there is some uncertainty in the capacitance 

arc using a 6900 bar transducer at a maximum pressure of of the system, which would alter the transducer calibration 

150 bar, and it is known that over such a wide pressure factor, ft would obviously be desirable to do a hcad-to-tail 
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in situ calibration of the transducer set-up, but neither 
time nor money to perform this has been available to date. 

The code also makes predictions of where the pump-tube 
piston ends up in the converging high-pressure section of 
the gun. The real piston initially has a conical hole in its 
forward end, as shown in figure 9(a). By the time the pis- 
ton ends up jammed in the high-pressure section of the 
gun, this hole is very much elongated and the forward end 
of the piston typically separated into several fingers. We 
have measured the final dimensions of two pistons and 
from these have determined the axial stations (in the 
co-ordinate system of fig. 8) of (1 ) the ends of the solid 
parts of the pistons and (2) the ends of the plastic fingers. 
Figure 14 is a sketch of one-half of the high-pressure sec- 
tion of the gun, in which are shown the experimentally 
measured stations of the piston data and the CFD pre- 
diction of the final position of the end for the one- 
dimensional solid piston. (The two experimentally mea- 
sured positions of the ends of the fingers fall on top of 
each other in fig. 14.) The code prediction is excellent, in 
that it falls roughly mid-way between the experimentally 
observed ends of the solid part of the piston and the ends 
of the fingers. 

The comparisons of the CFD and experimental data for 
the Ames 1 .5 in. gun presented up to this point have been 
for one operating condition, that of the average of shots 
600-607 in table 2. In the current test series, the gun has 
been fired at several other operating conditions (i.e., for 
shots 613-618) and data for these shots (excluding 
shot 617) are presented in table 2. Below, we make some 


limited comparisons from the data for the benchmark 
mean condition (“AV. 600-607”) and shots 613-618 in 
table 2. (We note here that for shots 616-618, two 
sections of pump tube were removed, reducing the pump 
tube volume by about 33 percent and the pump-tube fill 
pressure was increased to maintain the same mass of 
hydrogen as for shots 600—607.) Figure 15 shows the 
comparison of the CFD and experimental piston velocities 
for these six conditions, along with the ideal perfect 
agreement line. The agreement is very good. The 
predictions are, of course, tuned to agree at the 685 m/sec 
point. For the other five data points, the errors are 
0.8-1 .0 percent (6-7 m/sec) for three points and about 
1 .5 percent (1 1 m/sec) for two points. 

Figure 16 shows the comparisons of the CFD and experi- 
mental projectile muzzle velocities for the five of the six 
test conditions for which piston velocity data is shown in 
figure 15. The comparison for shot 615 is not shown, for 
the following reason. For this shot, for reasons of avail- 
ability, a diaphragm much thicker than usual was used. 
Upon disassembling the gun after the shot, it was apparent 
that, after rupture, the thicker diaphragm throttled the flow 
significantly just aft of the initial position of the projectile 
base. This is believed to be the reason for the very low 
experimental muzzle velocity for this shot (6.3 km/sec 
versus a CFD prediction of 6.95 km sec). Hence, the data 
for this shot is not shown in figure 16. For the usual 
diaphragms, there is no throttling effect, since adequate 
diaphragm relief is machined into the barrel for 1-2 cal- 
ibers downstream of the diaphragm station. 
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CFD projectile velocity, km/sec 

Figure 16. CFD and experimental projectile velocities for the Ames 1.5 in. gun. 


The point at 6.73 km/scc is the point for the condition 
“AV. 600-607” in table 2. Agreement between theory and 
experiment is almost perfect for this point, as has been 
previously noted. The two points with experimental veloc- 


ities of 6.82 km/sec (for shots 613 and 614) show velocity 
errors of about 3 percent (0.1 8 km/sec). For these same 
two points, the experimental piston velocities were about 
1 percent low (see fig. 13). The muzzle velocity is known 
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to he extremely sensitive to the piston velocity and the 
1 percent error in the latter velocity would lead to muzzle 
velocities 2-3 percent low. The remaining two data points 
show errors of 0. 1 percent and 1 .8 percent. Overall, the 
agreement between the CFD predictions and the experi- 
mental data was judged to be very good. 

5.4 Code Validation with Data from the Ames 0.28 in. 
Gun 

Data for a number of shots with the Ames 0.28 in. gun are 
given in table 4. A wider variety of gun operating 
conditions are represented here than for the data for the 

1.5 in. gun (see table 2). Piston masses of 155, 205, and 
222 grams and powder loads of 34, 40, and 50 grams were 
used. Three of the shots (shots no. 592, 594, and 596) had 
an extra diaphragm inserted in the pump tube to attempt to 
improve gun performance (see ref. 1 1). Projectile and pis- 
ton velocities were not successfully obtained for all shots. 
The first CFD modelling using the present code was done 
for shots 577-579. However, the code contained an error 
at this point which caused an error in the projectile veloc- 
ities, but did not affect the piston velocities. Hence, the 
comparison of projectile velocities is omitted for these 
three shots. The break valve rupture pressures shown in 
table 4 were calculated using the methods of reference 40. 
For the corresponding CFD calculations, these pressures 
were doubled in an attempt to allow for increases in yield 
and ultimate strength due to high strain rate effects (sec, 
for example ref. 4 1 ). This is now no longer our preferred 
method of calculating rupture pressures of the break 
valves. We now use, instead, the method of reference 36 
(as mentioned earlier) and use the value so calculated 
directly (without doubling) in the CFD calculations. 


Figure 17 shows the comparison of the experimental ver- 
sus CFD piston velocity data from table 4. Overall, the 
agreement, over a wide range of operating conditions, is 
very good. Eight data points are shown; for five of the 
points, the errors arc less than 2.5 percent, one point has 
an error of 3 percent and two points have errors of about 
4 percent. 

Six sets of projectile velocity data are shown in table 4. 
For two of these shots (shots 592 and 593), the experi- 
mental velocities are very low, nearly I km/sec below the 
CFD predictions. These two projectiles were known to 
have somewhat undersize diameters and it is believed that 
this resulted in massive gas blow-by past the projectiles 
and, hence, the low muzzle velocities. The projectile 
diameters for shots 584 and 594-596 were much closer 
fits to the barrel diameter. The projectile velocity data for 
these four shots are plotted in figure 18. The experimental 
values run, in general, about 4 percent below the CFD 
values. To allow for this, we have plotted the line 
(experimental muzzle velocity) = 0.958 x (CFD muzzle 
velocity) on the graph (dashed line). The data points fall 
within 2 percent of this line, indicating that the code pre- 
dicts the variations of gun performance quite well, once 
the overall 4 percent difference is allowed for. The 
0.28 in. gun has a history of not shooting as well as other, 
larger, guns at Ames. This may, perhaps, be due to the 
occurrence of a certain amount of gas blow-by even for 
the best projectiles. This, is turn, may be due to the 
inevitably poorer machining tolerances and finishes which 
will occur, relative to the gun size, for a small gun, since 
these tolerances and finishes lend to have fixed absolute 
values. 


Table 4 Experimental and CFD data for shots on NASA Ames 0.28 in. gun 


Shot 

no. 

Powder 

mass 

(gm) 

Piston 

mass 

(gm) 

Proj. 

mass 

(gm) 

Break 

valve 

pressure 

(bar) 

Hydrogen 

pressure 

(bar) 

Pump 
tube 
dia p. 
(bar) 

Pump 

tube 

length 

m 

CFD 

piston 

velocity (m/s) 

Expcr. 

piston 

velocity (m/s) 

CFD 

project. 

velocity 

(km/s) 

Expcr 

project. 

velocity 

(km/s) 

577 

34 

204 

0.334 

517 

1.52 


100 

665.8 

666 



578 

40 

205 

0329 

517 

1.52 


100 

729.7 

748.2 



579 

40 

155 

0.324 

517 

1.52 


100 

819.5 

800.7 



583 

50 

221 

0.391 

690 

1.66 


100 

790.7 

822.3 



584 

50 

222 

0.371 

690 

1.66 


100 



7.134 

6.741 

585 

50 

222 

0.398 

690 

1.66 


100 

785 

822.1 



592 

40 

206 

0.339 

690 

6.9 

58.3 

75 

684.5 


6.65 

5.623 

593 

40 

204 

0.335 

690 

2.16 


75 

710.5 

706.5 

6.553 

5.73 

594 

40 

206 

0.323 

690 

6.9 

28.3 

75 

697.6 


6.885 

6.714 

595 

40 

207 

0.336 

690 

2.16 


75 

710.5 

688.5 

6.553 

6.346 

596 

40 

206 

0.336 

690 

3.92/1.01 

28.3 

75 

700.5 

710.1 

6926 

6.504 
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5.5 Summary of Code Validation with Actual Gun 
Data 

In this section, we sum up the validation of the code 
against experimental data from the Ames 1.5 in. and 
0.28 in. guns. The validations were presented in detail in 
the two preceding sections. The first set of validations 
were for data at the benchmark operating condition of the 
Ames 1.5 in. gun (shots 600-607 in table 2). The valida- 
tions consisted of comparisons of the following CFD and 
experimental data: 

(1) pressure histories in the powder chamber, 

(2) pressure histories in the pump tube, 

(3) projectile muzzle velocity, and 

(4) final position of end of pump-tube piston. 

As described, agreement between the CFD and experi- 
mental data was excellent for all of these variables. 

The second set of validations consisted of comparing CFD 
and experimental piston and projectile velocities for a 
number of different piston masses and powder loads for 
two very different guns — the Ames 1 .5 in. and 0.28 in. 
guns. Data for a greater range of operating conditions 
were available for comparisons with the 0.28 in. gun. 

Very good agreement was found between theory and 
experiment for piston velocities. Projectile velocities were 
very well predicted (errors of 1-3 percent) for the Ames 

1.5 in. gun. For the Ames 0.28 in. gun, the experimental 
projectile velocities were about 4 percent low. The small 
gun has a history of relatively poor performance and it is 
felt that this may be due to gun tolerance problems and 
possible blow-by of gas past the projectile which is not 
modelled by the code. 

Overall, the validation of the code against actual gun data 
was judged to be very good. 

5.6 Other Predictions of CFD Code 

In this section, we present a number of code predictions 
which illustrate the usefulness of the code in (1) under- 
standing phenomena within a two-stage gun and 
(2) choosing more favorable gun operating conditions (for 
example, lowering maximum gun and projectile base 
pressures while maintaining muzzle velocity). We are not 
currently able to make any direct measurements within 
certain key regions of our guns (i.e., high-pressure section 
and barrel), so the code’s ability to give an “X-ray” pic- 
ture of the internal fluid dynamics of the gun is very valu- 
able for optimizing gun operating conditions. Ail of the 
predictions discussed in this section are from the CFD 
calculation for the benchmark gun operating condition of 
shots 600-607 in table 2. 


Figure 19 shows CFD position and velocity histories of 
the front of the pump-tube piston. The piston accelerates 
in an s-shaped curve up to its maximum velocity of about 
700 m/sec and then decelerates much more rapidly in the 
high-pressure section of the gun. As the piston enters the 
contraction section of the gun and is decelerating as a 
whole, there is some tendency for re-acceleration of the 
front end of the piston, due to the rapid area reduction. 

For the case of figure 19, this appears only as a small kink 
at a velocity of 540 m/sec, but for other, more violent 
operating conditions, this can result in a significant 
re -acceleration, for example, from 700 to 1000 m/s ec. 

Figure 20 shows the CFD pressure histories at two sta- 
tions (3355.033 and 3380.127 cm) in the high-pressure 
section of the gun. These stations are located at about 
35 percent and 65 percent of the way along the contrac- 
tion (see fig. 8). A series of waves reflects between the 
front of the piston and the diaphragm (or the projectile 
base, after diaphragm rupture). These can be seen in the 
pressure histories of figure 20, up to pressure levels of 
about 2000 bar (also in figs. 1 2 and 13). Note the very 
high (6000-7000 bar) pressures calculated to occur in the 
high-pressure section of the gun. Obviously, limiting the 
maximum pressures reached here (and elsewhere in the 
gun) would be very desirable to obtain long component 
lifetime. By a judicious selection of gun operating condi- 
tions, it is possible to reduce the maximum pressures, 
while maintaining muzzle velocity, or, conversely, muzzle 
velocity can be increased while maintaining the maximum 
pressure levels in the gun. 

Figures 21 and 22 show pressure snapshots at five differ- 
ent times during the gun-firing cycle. The two figures are 
identical except for the expanded time scale of figure 22. 
The three vertical dashed lines represent the position of 
the contraction of the high-pressure section of the gun. 

The outside two lines represent the beginning and end of 
the contraction (as modelled in the CFD code, see 
sec. 5.1) and the middle line in the slight break in the 
slope of the contraction. The intervals between the snap- 
shots are not equal, but are about 0.6, 0.3, 0.3, and 
0.4 msec. The projectile base is at the right end of the 
curves. We note again the very high maximum pressures 
(6000-7000 bar) for the second and third curves (c.f., with 
fig. 20). For the first four curves, there is a roughly linear 
slope on the left-hand side of the curves. This represents 
the piston, with the rear of the piston being at essentially 
zero pressure and the front face of the piston at the maxi- 
mum pressure point. The roughly uniform pressure gra- 
dient over the length of the piston causes the piston to 
slow down very rapidly (see fig. 19). Note that the snap- 
shots in the time period of the highest pressures show a 
pressure drop by about a factor of two from the high- 
pressure reservoir to the start of the barrel (at about 
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Figure 20. CFD predictions of pressure histories in high-pressure section of gun. 


station 3410 cm). This is due to the start of the barrel act- 
ing as a sonic throat. The pressure at the projectile base is 
calculated to be much less than the maximum pressures in 
the high-pressure section of the gun (i.e., 1000-1500 bar 
maximum versus 6000-7000 bar). Finally, in the first four 
curves and, in particular, for the third and fourth curves, 
pressure waves can be seen in the barrel running between 
the piston front and the projectile base. 


Figure 23 shows the pressure histories in the most forward 
hydrogen cell (solid line) and at the projectile base (dotted 
line). This hydrogen cell in question is adjacent either to 
the break valve (diaphragm) of the projectile base. The 
two curves fall on top of each other for much of the lime 
interval shown. A series of waves of increasing ampli- 
tudes oscillates back and forth between the front of the 
piston and the break valve until the diaphragm breaks at 
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Figure 21. CFD predictions of pressure snapshots at several times in the gun. 
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Figure 22. CFD predictions of pressure snapshots at several times in the gun. (As for fig. 21, but with expanded time 
scale.) 


about 57.7 msec. Thereafter, the waves continue to oscil- spikes are separated by regions of much lower pressures 

late between the front of the piston and the base of the (300-400 bar). This can be shown to be due the shock 

projectile until about 60.1 msec, after which the projectile focusing action of the contraction section of the gun, 

base pressure falls continuously. The effective average which has an average full angle of about 8.7 deg. Using a 

projectile base pressure over the length of the barrel is much larger angle, such as 40 deg or 60 deg, will make 

only about 700 bar, although with many spikes, the high- the pressure history at the projectile base much more 

est of which reach about 1600 bar. These very sharp uniform. 
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Figure 23. CFD predictions of pressure history at projectile base and in the last (most forward) hydrogen cell (cell 3,53). 
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Figure 24. CFD predictions of projectile base pressure as a function of projectile position in barrel. 


Figure 24 shows the projectile base pressure plotted as a can he readily be determined as the mean pressure in lig- 
lunction of the position of the projectile in the barrel. The urc 24, over the length of the barrel, 

same waves shown in figure 23 appear also in figure 24. Figurc 25 shows CFD fK)sitjon and vc | ocily histories of 

From figure 24, one may determine the position of the the pro j ec tile. The steepest slopes of the velocity curves 

projectile in the barrel when .1 was struck by the various corrcspond lo thc arrival, at the projectile base, of the 

pressure waves. Also, since (work) = (pressure) x (area) x surc ks shown in ngurcs 2 3 and 24. 

(distance), the effective average projectile base pressure 
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Figure 25. CFD predictions of velocity and position histories of projectile . 


5.7 Discussion of Anomalous Piston Behavior 

In this section, we discuss what appears to be very 
anomalous behavior of the pump tube piston. The piston 
shank diameter is considerably smaller than the tube bore 
(see figs. 8 and 9) being designed to ride free of the bore 
to reduce piston friction. However, following the model of 
section 2.5 and Appendix C (based on that of ref. 4), the 
shank should bear against the tube bore for part of the 
piston stroke cycle for the following reasons. The maxi- 
mum powder-burn pressure is approximately 1000 bar 
(fig. 1 1) and portions of the shank should see axial pres- 
sures of up to 700 bar. The maximum hydrogen pressures 
are about 7000 bar (figs. 20-22) and portions of the shank 
should see axial pressures up to 4000-5000 bar. Now, the 
low strain rate yield strength and tensile modulus of elas- 
ticity (E) of polyethylene are given in a number of refer- 
ences (e.g., ref. 42) as roughly 200 bar and 9000 bar, 
respectively. From the longitudinal and shear sound speed 
and density data of reference 21 for polyethylene, using 
the equations of reference 43, we can deduce an E value 
of about 27,000 bar and a Poisson’s ratio of about 0.4. 
Now, since the axial pressures in the shank very much 
exceed the stated yield strength at low-strain rate, based 
on the latter numbers, at least part of the shank should 
collapse and bear on the tube wall. Even if we assume 
that, at high strain rate, the yield strength of the piston is 
somehow very much increased above the low strain rate 
values, simply based upon elastic expansion of the piston 
under the applied axial pressure, the piston should expand 
laterally and bear against the tube wall. 


However, upon examining a number of pistons after fir- 
ing, we observe the following. The rear land clearly rides 
along the bore (as it must), and evidence of this is as 
follows. 

(1) The land diameter is within 0.002-0.007 cm of 
that of the bore (some measurements being very slightly 
larger and others, very slightly smaller than the bore). 

(2) The original circumferential machining marks 
are heavily worn down and, in some cases, removed. 

(3) The land shows scratches and streaking in the 
axial (travel) direction. 

On the other hand, the shank (except where it has been 
jammed into the contraction of the high-pressure section) 
shows absolutely no evidence of riding on the bore as 
indicated by the following observations. 

( 1 ) The shank diameter is consistently about 0.20 (or 
more) cm less than the bore diameter, 

(2) The original circumferential machining marks 
remain undamaged. 

(3) The shank shows no sign of axial scratches and 
streaking (except at the whisker gauge azimuth). 

Thus, the piston appears to act as though both of the fol- 
lowing are true: 

(I) Its high strain rate yield stress is much higher 
than the quoted low-strain rate value. 
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(2) Its high strain rate tensile/compressive modulus 
is much higher than the quoted low-strain rate value (or 
alternatively that the high strain rate Poisson’s ratio is 
much lower than the value of 0.4 derived from the data of 
ref. 21). 

As noted in section 3.5, the piston would be calculated to 
bounce back after maximum gas compression if the be at 
the front of the piston were not changed to stop the piston 
motion when the velocity of the front of the piston reaches 
zero. This was done in order to agree with the observed 
experimental behavior. This jamming behavior suggests 
that the high strain rate shear yield stress may be substan- 
tially higher than that estimated from the quoted low- 
strain rate tensile yield stress as described in Appendix C. 
If the high strain rate shear yield stress (and the friction 
coefficient) are sufficiently high, this would suffice to halt 
the piston naturally, without having to change the be. The 
observed jamming behavior of the piston is thus, in a gen- 
eral way, consistent with the anomalous yielding and 
deflection behavior discussed above for the shank of the 
piston. 

From the above discussion, there appear to be some very 
serious unresolved issues regarding high strain rate 
behavior of the pump tube piston and the modelling of 
same. We hope that the above discussion may perhaps, to 
a small degree, aid and stimulate the ballistics community 
towards the development of better modelling of piston 
phenomena. 

6 Summary and Conclusions 

We have presented a new code for the calculation of the 
performance of two-stage light gas guns. This code is 
based on the Godunov method and is third-order accurate 
in space and second-order accurate in time. The Riemann 
solver used is exact for shocks and uses a very accurate 
power law integration for expansion waves. Realistic EOS 
are used for all media. The code includes modelling of 
friction and heat transfer for powder gas, hydrogen, the 
pump-tube piston and the projectile. A simple nonequilib- 
rium turbulence model is included for the gas flows and 
the predictions of skin friction and heat transfer to the 
tube walls in the gas flows are modified accordingly. 
Gunpowder bum in the first-stage breech is modelled 
using standard ballistic techniques. 

The code was first validated with a number of analytical 
solutions. These included (1) Riemann’ s shock tube 
problem at pressure ratios up to 10 6 , (2) impact of two 
polyethylene plates at a closing velocity of 20 km/sec, 

(3) flow through a convergent-divergent supersonic noz- 
zle with an area ratio of 16 to 1, and (4) gunpowder burn 


in a closed bomb. All these code validations were found 
to be very good to excellent. 

The code was then validated by comparing its predictions 
with experimental data from the Ames 1.5 in. and 0.28 in. 
light gas guns. These data included: 

( 1 ) powder-chamber pressure histories, 

(2) pump-tube pressure histories, 

(3) piston velocities, 

(4) projectile velocities, and 

(5) the observed final position of the front of the 
pump-tube piston. 

To model the powder-chamber pressure histories and the 
piston velocity properly, it was found to be necessary to 
do the following: 

(1) increase the burning rate over that given by the 
powder manufacturer by 60-90 percent and 

(2) reduce the pump-tube piston friction below that 
predicted by direct, unadjusted application of the piston 
friction model. 

Reasons for these changes and the way in which the 
parameters were modified were discussed in depth. After 
the powder burn rate and the piston friction were adjusted 
for one standard operating condition of each gun, they 
were then held constant for all other operating conditions. 
Predictions of the CFD code and the experimental mea- 
surements were in very good agreement for two very dif- 
ferent guns (the Ames 1 .5 in. and 0.28 in. guns), over a 
considerable range of operating conditions (particularly 
for the 0.28 in. gun). 

Predictions of the code for the piston and projectile posi- 
tion and velocity histories for the 1 .5 in. gun were pre- 
sented. Calculated pressure histories and snapshots in the 
high-pressure section of the gun and at the projectile base 
were also presented. Discussion of these histories showed 
their usefulness in optimizing gun operating conditions, 
for example, to maximize muzzle velocity while maintain- 
ing given maximum gun pressures or, conversely, to 
reduce the maximum gun pressures while maintaining 
muzzle velocity. 

Finally, observations made on pump tube pistons after fir- 
ing strongly suggest that the high strain rate yielding and 
deflection behavior of the pistons in the pump tube is very 
different than that which would be predicted using low- 
strain rate data. It is suggested that considerably more 
work in this area may be necessary to allow improved 
predictions of piston behavior to be made. 
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APPENDIX A 


Gas-Phase Friction and Heat Transfer 

We start this analysis with well known (ref. 18 and 19) 
skin friction correlations for pipe flow. To correct for the 
effects of Mach number (M) and the difference between 
the wall temperature (T w ) and the average stream temper- 
ature (T), we use the correlations developed by van Driest 
(ref. 20). We have developed a “reference temperature” 
technique which reproduces the variations of skin friction 
coefficient, f (or Cf), with M and T w /T given by van 
Driesfs correlations within about 10 percent over a wide 
range of parameters without requiring the iterative solu- 
tion of equations (66) or (71) in reference 20. We empha- 
size that our “reference temperature” is not the same as 
the usually used reference temperature (ref. 44). Our 
“reference temperature” (T) is given by 


T = 0.9T + 0.03M 2 T + 0.46T W (A-l) 


We then evaluate the density and viscosity at the reference 
conditions 


P' = P^ (A-2) 

p , = fi , (T\p < ) (A-3) 


where p is the mean cell density. The basic, low pressure 
viscosities as functions of temperatures are fit using sim- 
ple power laws to data from Golubev (ref. 45). This refer- 
ence gives data for H2, N2, CO, CO2, and H2O. For the 
hydrogen pump-tube gas, the H2 data is used directly. For 
the powder gas, the composition of the powder gas was 
taken from data sheets (ref. 46) provided by the powder 
manufacturer. The viscosity of the powder gas mixture at 
various temperatures was then found by adding up the 
individual viscosities from reference 45 multiplied by the 
mole fractions. The resulting gas mixture data was then fit 
with a power law curve. Golubev (ref. 47) gives density 
corrections for viscosity for H2, N2, and CO2 These cor- 
rections as a function of density were fit with a sum of 
two terms with density ratios to two different powers and 
used to estimate a density correction for hydrogen and 
powder gas. Lacking pressure correction data for CO and 
H2O, the data for N2 and CO2 was assumed to apply to 
CO and H2O, respectively. 


With p‘ and |T determined, the Reynolds number (Re) for 
the cell in question is calculated as 


Re ' d 


p 1 uD 
V' 


(A-4) 


where u is the cell-center gas velocity and D is the mean 
diameter of the cell. The skin friction coefficient (ref. 18) 
for pipe flow as a function of Re is then fit with the fol- 
lowing expressions 


r = 0.049(Re' D ) °- 2 , 

Re' D > 5507 

(A-5) 

r = 0.00875, 

5507 > Re' D > 1828 

(A-6) 

r- 16 

1 828 > Re’ D 

(A-7) 

Re’ n 


for the turbulent, transition and laminar regimes, respec- 
tively. The skin friction coefficient is then corrected for 
nonequilibrium turbulence effects as follows 


f 

f TC = r LAM + ( r _r LAM ) 

V 


TKEncq 

tke~ 


,0.5 


(A-8) 


where TKE e q and TKEpeq are the equilibrium (steady- 
state) and nonequilbrium turbulent kinetic energies for the 
cell, Flam is the laminar skin friction coefficient, always 
calculated from equation (A-7) and fjc the cor ~ 
rected skin friction coefficient. TKE neq is calculated as 
shown in Appendix B. fjc is referenced to the density at 
the reference temperature; it is modified to refer to the 
mean cell density as follows 


f = r 



(A-9) 


The wall-friction force on the gas, Fjy, for the cell in ques- 
tion is then given by 


Ffr = — jfpu|u|nDAx (A- 10) 
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where the wall area for the cell is tcDAx, with Ax being 
the cell length. Absolute value signs are required on one 
u term to maintain the proper sign of Ff r . To calculate the 
heat transfer from the gas to the walls, Reynolds’ analogy 
(ref. 21) is used in the following form, which is correct if 
the Prandtl number is taken to be unity: 


Combining equations (A* 1 1)— (A- 1 3) and using the rela- 
tion between h, e, p, and p, we obtain, finally, 



(A- 14) 


Q w =^(H-h w ) (A- 1 1) 

u 

where H is the mean cell total enthalpy, h w is the static 
enthalpy evaluated at the wall and Q w is the heat-flow 
rate from the cell to the wall. The mean cell total enthalpy 
is related to the mean cell static enthalpy, h, by the usual 
relation 


H = h + — (A- 12) 

2 

Since h w is not directly available as the cell center compu- 
tational fluid dynamics (CFD) solution progresses, it is 
approximated as follows 


which is used to calculate the heat-flow rate from the cell 
to the wall in gas-flow regions. 

For regions with gas and solids flowing together 
(gunpowder/powder gas), Ff r and Q w are based on the 
density, enthalpy and sound speed of the gas only. The 
presence of solids is ignored. This two-phase flow situa- 
tion applies only in the upstream part of the pump tube for 
the first part of the solution before all the powder is 
burned. Currently, T w is prescribed and held fixed at 
300 K. Estimates of the wall heating at various locations 
in the gun have been made. In most of the pump tube and 
the downstream part of the barrel, these effects are rather 
small. However, as is well known, in the high pressure 
coupling and in the upstream part of the barrel, wall heat- 
ing is significant. However, the results of the code 
obtained to date are sufficiently good and useful even 
without wall-heating effects to warrant presentation in the 
current publication. At a later time, it is intended to mod- 
ify the code to allow for wall-heating effects. 



(A- 1 3) 
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APPENDIX B 

Gas-Phase Nonequilibrium Turbulence Model 

A simple model is developed which assumes that the 
nonequilibrium turbulence kinetic energy (TKE) relaxes 
towards the equilibrium value for the flow in question 
(TKE eq ) with an e-folding length (L^) which is a certain 
number of tube diameters. (The e-folding length is the 
length over which the difference between the nonequilib- 
rium and equilibrium TKE will relax to 1/eth of its origi- 
nal value in a steady, constant area flow.) Hinze (ref. 22) 
presents an extensive discussion of the fully developed 
low speed turbulent pipe flow measurements of Laufer 
(ref. 23). We estimate L e using (1 ) Laufer’s graphs 
(presented in Hinze) of the TKE distribution across the 
pipe radius, (2) Laufer’s graphs (also presented in Hinze) 
of the TKE production distribution across the pipe radius, 
and (3) Schlichting’s (ref. 48) values for the ratio of max- 
imum to mean velocity for low speed, fully developed 
pipe flow. From these data for Reo = 5 x JO 5 , we estimate 
L c = 3.27 x (pipe diameter). The range of Re for hydrogen 
flow in the pump tube is typically 3 x 1(P to 3 x 10 7 . The 
Re for the data of references 22 and 23 is within our 
range, but towards the low end of it. However, turbulent 
pipe flow docs not appear to change very rapidly with Re 
over the Re range of interest (at least over the range 
3 x 10 5 to 3 x 10 6 reported in ref. 24). Hence, we use the 
value of L e given above as a rough estimate in our CFD 
model. The relaxation term in our model thus becomes 


d(TKE) = — (TKE eq -TKE) (B-l) 


where d(TKE) is the change in TKE which takes place 
when the flow moves a distance dx, and we use 

TKE eq =0.00929u 2 (B-2) 

also taken from the data of Laufer (ref. 23) for 

Rep = 5x10^. For simplicity in the equations, we have 

dropped the subscript “neq” from TKE nC q’, i.e,, “TKE” in 


the present equations corresponds to TKE ne q in 
Appendix A. For one timestep dt, the distance that the 
flow moves is simply udt. Since the tube changes diame- 
ter in the gun model, L e is not fixed, but is taken to be 
equal to RlD, where Rl = 3.27, as discussed above. 
Inserting these two results into equation (B-I) yields the 
following equation for the relaxation term of the TKE 
equation 

d(TKE) = (TKE eq - TKE ) (B-3) 


A difficulty with equation (B-3) is that, as it stands, there 
will be no TKE relaxation if the velocity, u goes to zero. 
Since the TKE will obviously relax due to the turbulent 
motion itself, even if u = 0, we have modified equa- 
tion (B-3) by replacing u with (u 2 + 2TKE) 0 * 5 to yield our 
final form for the relaxation term of the TKE equation, as 
follows 


r~2 

d(TKE)= — J r 2TKE dt(TKE eq - TKE) (B-4) 


To calculate the changes of TKE within any cell over a 
timestep, equation (B-4) is used, along with the usual 
terms taking account of convection of TKE across the cell 
boundaries of cells. The nonequilibrium value of TKE 
thus calculated is then used to modify the skin friction 
coefficient as described in Appendix A. 

Examination of the nonequilibrium and equilibrium TKE 
values calculated by the code shows the following. For the 
acceleration of the piston and the projectile, the TKE val- 
ues are calculated to lag only slightly behind the equilib- 
rium values, even during the very high initial accelera- 
tions of the piston and projectile. When the piston slows 
down very rapidly in the high-pressure section of the gun, 
however, the TKE in the rapidly decelerating gas in front 
of the piston is calculated to rise well above the equilib- 
rium values. The inclusion of the nonequilibrium TKE 
model was found to make only minor changes 
( 1-3 percent) in the code predictions of pressures and 
piston and projectile velocities. 
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APPENDIX C 


p x = axial pressure (from CFD code results) 
s n = stress in solid normal to the x direction 

Solid-Phase Friction and Heat Transfer 


The method used is a variation of that presented in 
reference 4. We begin by calculating the normal 
stress between the solid body (piston or projectile) 
and the tube wall (a n ) using the following logic and 
equations. Equations (C-2HC-7) are written in the 
form of two nested block IF FORTRAN statements. 




(C-I) 


IF4>>H 


°n = Px _ °y 

(C-2) 

ELSE 


< 

II 

1 

(C-3) 

IF <t» < |i 


HP X 

Or, = Gy + 

n 1 l-n 

(C-4) 

ELSE 



f \ 


-E 


+ W* 


_ _ v ip ; 


(C-5) 

n " l-F 


ENDIF 



If(Px ” ) — ^y * 

“ Px ^y 

(C-6) 

“ Px ) — ^y * 

= Px + ^y 

(C-7) 


In the above equations, since p x is positive for pressures 
in the CFD code, we have taken a n to be positive for 
compression, which is the reverse of the usual strength of 
materials convention. 

The above set of equations looks formidable but actually 
embodies a rather simple elastic-plastic model as follows. 
First, the normal stress on the solid (a n ) is calculated 
assuming elastic behavior, taking into account the applied 
pressure in the x direction (p x ) and the initial jamming of 
the solid rod into the tube. If the difference lo n - p x l so 
calculated is less than the yield stress, C y , this value of a n 
is used unmodified. If I a n - p x l > o y , then a n is set equal 
to p x ± Oy , as appropriate. This represents the plastic 
conditions part of the model. The case where the solid rod 
is smaller than the tube (free), either initially (unstressed) 
or after the application of p x , is allowed for in the above 
logic set. 

Again following reference 4, the normal stress is limited 
as follows (so that it cannot become negative): 


o n = max(0,a n ) 


(C-8) 


A number of references 25-28 have discussed high speed 
models for dynamic sliding friction. Reference 26 points 
out that the experimental data of reference 25 for nylon on 
steel up to about 0.7 km/sec can be fairly well fit with a 
curve of the form 


(i d = Au 


-0.4 


(C-9) 


ENDIF 

where: 

r t = tube radius 

Tp = unstressed initial radius of solid material 
(outside tube) 

G y = yield stress of solid 

E = Young’s modulus of solid 

|i = Poisson’s ratio of solid 


where u is the relative velocity, A is a constant and p d is 
the dynamic friction coefficient. Reference 26 used equa- 
tion (C-9) for studies of model wear in a two-stage gas 
gun. Reference 27 used the same equation for friction 
modelling in an electromagnetic launcher. In reference 28, 
equation (C-9) for the friction coefficient was generalized 
to the following expression for pf, at any velocity 
(including zero velocity), 

Pf = min|p s , Au -h j (C-10) 


PRECEDING PAG^ M _NO S; F,LMED iv 4, 



where |i s is the static friction coefficient. We have used 
equation (C-10) recommended by reference 28 and have 
reevaluated the coefficients A and b directly from the data 
of reference 25 for nylon on steel. We have obtained 
b = 0.4224 and A = 8.377 (if the velocity is in cm/sec). 
From figure 242 of reference 25, |i s is estimated as 0.185. 
Lacking dynamic friction data for piston and projectile 
materials such as polyethylene and Lexan, we have simply 
scaled the A coefficient in equation (C-10), obtained from 
the data of reference 25 for nylon on steel, by the 
respective static friction coefficients to make first esti- ' 
mates for the dynamic friction coefficients for the other 
materials. 

With an estimate for the friction coefficient available, the 
wall-shear stress, x w , is calculated from the following two 
equations 


*w — M-f^n 


(C-ll) 


t 


w 


min 





(C-12) 


stress to the shear yield stress, here taken to be 1 / V3 
times the tensile yield stress, following reference 29. With 
x w available, the wall friction force on the solid in a given 
cell can be found simply by multiplying x w by the area of 
the cell in contact with the wall. The heat- flux rate to the 
solid, q w , due to the solid friction work is simply taken to 
be one half the shear stress times the velocity, as follows 


— _ XijuU 


(C-13) 


This term, which appears as a source term in the energy 
equation, represents energy loss from the solid to the wall 
due to frictional heating occurring at the wall. The other 
half of the heat generated by the frictional work is 
assumed to flow to the solid and does not represent an 
energy loss from the cell and therefore should not be 
included in the energy equation. If a point mass projectile 
is assumed, the wall shear stress and heat-flux source 
terms are calculated midway along the corresponding 
(equal mass) real length projectile and multiplied by the 
proper wall-projectile contact area for the real length 
projectile. 


Equation (C-l 1) is simply the basic friction relation and 
equation (C-12) limits the maximum possible wall shear 
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